Carotenoids and immune data

Introduction

Statistical analysis of plasma carotenoids, plasma cytokines and immune cells from randomized cross-over USDA inflammation clinical trial. Subjects consumed both low lycopene tomato (yellow) and high lycopene tomato-soy juices (red) for 4 weeks each.

Crossover clinical trial design supplementing individuals with obesity 360 mL of a low carotenoid tomato juice or a high lycopene tomato-soy juice daily. Daily serving of low carotenoid tomato juice consisted of ~1.5mg lycopene/day while high lycopene tomato-soy juice intervention consisted of 54 mg lycopene/day in addition to 210 mg total soy isoflavones/day.

Load libraries

library(tidyverse) # data wrangling
library(readxl) # read in excel files
library(janitor) # clean up names in dataset
library(corrr) # finding correlations
library(rstatix) # stats
library(lme4) # mixed linear modeling
library(knitr) # aesthetic table viewing
library(lmerTest) # add pvalue column to lmer models
library(purrr) # create functions
library(broom.mixed) # generate tidy data frames for lmer results
library(MuMIn) # lmer model testing using AICc
library(kableExtra)
library(ggthemes)
library(ggtext)
library(ggpubr)

Read in data

# load data
meta_table <- read_excel("CompiledData_Results_Meta.xlsx",
                         sheet = "metadata_corrected_withsequence")

# clean up variable names 
meta_table <- clean_names(meta_table)

str(meta_table)

## tibble [60 × 89] (S3: tbl_df/tbl/data.frame)
##  $ patient_id                          : num [1:60] 6101 6102 6103 6104 6105 ...
##  $ sex                                 : chr [1:60] "F" "M" "F" "M" ...
##  $ age_at_enrollment                   : num [1:60] 58 65 60 54 40 57 68 75 36 46 ...
##  $ bmi_at_enrollment                   : num [1:60] 31.1 36.8 30.3 33.1 30.4 ...
##  $ intervention                        : chr [1:60] "Red" "Red" "Red" "Yellow" ...
##  $ sequence                            : chr [1:60] "R_Y" "R_Y" "R_Y" "Y_R" ...
##  $ intervention_week                   : num [1:60] 2 2 2 2 2 2 2 2 2 2 ...
##  $ pre_post                            : chr [1:60] "pre" "pre" "pre" "pre" ...
##  $ period                              : chr [1:60] "B1" "B1" "B1" "B1" ...
##  $ if_ng                               : num [1:60] 1.16 11.74 0.93 0.63 1.88 ...
##  $ il_1b                               : num [1:60] 15.26 28.75 8.62 11.04 22.5 ...
##  $ il_2                                : chr [1:60] "0.98" "9.32" "0.55000000000000004" "0.65" ...
##  $ il_6                                : num [1:60] 2.18 5.83 2.53 1.29 2.25 ...
##  $ il_8                                : num [1:60] 2.86 4.35 3.93 1.93 4.14 3.05 1.35 2.01 4.43 2.86 ...
##  $ il_10                               : chr [1:60] "0.05" "0.71" "1.89" "7.17" ...
##  $ il_12p70                            : num [1:60] 2.86 60.93 5.46 2.79 4.06 ...
##  $ mcp_1                               : num [1:60] 51.9 189.3 224.9 244.7 167.9 ...
##  $ tn_fa                               : num [1:60] 14.9 92.3 26.5 19.6 49.6 ...
##  $ il_13                               : chr [1:60] "9.64" "92.04" "21.61" "23.17" ...
##  $ il_5                                : num [1:60] 3 13.04 2.74 5.45 5.38 ...
##  $ il_1ra                              : num [1:60] 3.52 11.36 3.71 3.38 6.08 ...
##  $ il_12p40                            : num [1:60] 51.5 139.4 81 78.7 175.9 ...
##  $ gm_csf                              : num [1:60] 10.9 195.7 17.9 27.4 72.4 ...
##  $ il_4                                : chr [1:60] "NA" "8.8000000000000007" "1.6" "4.49" ...
##  $ x01_cd45_cd66b_lymph_dc_mono        : num [1:60] 99.3 95.6 99.7 99.7 96.4 ...
##  $ x02_cd45_cd66b_grans                : num [1:60] 0.682 0.33 0.305 0.285 0.696 ...
##  $ cd45_cd14                           : num [1:60] 46.9 71 69.8 55.8 42.4 ...
##  $ cd45_cd20                           : num [1:60] 44.2 59.9 60.6 54 31.9 ...
##  $ x03_cd3_cd45_cd3_t_cells            : num [1:60] 25.7 43.7 54.5 14.2 14.3 ...
##  $ x04_tc_rgd_cd3_ab_t_cells           : num [1:60] 22 42.6 53.1 11.2 13.6 ...
##  $ x05_cd4_cd8_cd8_t_cells             : num [1:60] 11.96 3.84 13.23 6.5 5.72 ...
##  $ ccr7_cd8                            : num [1:60] 4.385 1.54 4.186 0.962 1.882 ...
##  $ ccr7_cd8_2                          : num [1:60] 7.56 2.28 9.04 5.52 3.83 ...
##  $ x06_cd45ro_cd45ra_naive_cd8         : num [1:60] 3.692 0.233 3.753 0.652 1.667 ...
##  $ x07_cd46ro_cd45ra_cm_cd8            : num [1:60] 0.506 0.955 0.263 0.196 0.146 ...
##  $ x08_cd45ro_cd45ra_em_cd8            : num [1:60] 6.555 0.395 7.505 4.528 1.807 ...
##  $ x09_cd45r0_cd45ra_te_cd8            : num [1:60] 0.826 1.322 1.2 0.775 1.661 ...
##  $ x10_cd38_hladr_activated_cd8        : num [1:60] 0.2543 0.0823 0.1523 0.1347 0.0742 ...
##  $ x11_cd4_cd8_cd4_t_cells             : num [1:60] 9.14 37.12 37.51 3.64 6.98 ...
##  $ ccr7_cd4                            : num [1:60] 7.98 33.92 32.87 2.49 5.3 ...
##  $ ccr7_cd4_2                          : num [1:60] 1.16 3.19 4.64 1.15 1.68 ...
##  $ x12_cd45ro_cd45ra_naive_cd4         : num [1:60] 3.163 3.77 17.23 0.693 1.731 ...
##  $ x13_cd45ro_cd45ra_cm_cd4            : num [1:60] 3.51 22.39 8.5 1.25 2.63 ...
##  $ x14_cd45ro_cd45ra                   : num [1:60] 1.053 2.964 3.919 0.893 1.488 ...
##  $ x15_cd45ro_cd45ra_te_cd4            : num [1:60] 0.0257 0.0664 0.5457 0.1154 0.092 ...
##  $ x16_cd38_hladr_activated_cd4        : num [1:60] 0.321 0.698 0.384 0.123 0.14 ...
##  $ cd45_cxcr5_th                       : num [1:60] 8.291 28.987 0.133 2.95 6.456 ...
##  $ x17_cd25_cd127_tregs                : num [1:60] 0.348 0.545 0.438 0.131 0.193 ...
##  $ x18_ccr4_cd4_total_ccr4_treg        : num [1:60] 0.317 0.494 0.316 0.123 0.18 ...
##  $ x19_cd45ra_cd45ro_ccr4_treg_naive   : num [1:60] 0.00973 0.03494 0.00676 0.00955 0.00744 ...
##  $ x21_cd45ra_cd45ro_ccr4_treg_memory  : num [1:60] 0.3017 0.4366 0.2905 0.0949 0.166 ...
##  $ x20_hladr_total_ccr4_treg_activated : num [1:60] 0.2867 0.2547 0.1017 0.0838 0.1119 ...
##  $ x22_cxcr3_ccr6_th1                  : num [1:60] 0.7981 2.0038 0.0209 0.4805 0.8362 ...
##  $ x23_cxcr3_ccr6_th2                  : num [1:60] 5.6202 13.4081 0.0747 1.2236 4.3881 ...
##  $ x24_cxcr3_ccr6_th17                 : num [1:60] 1.1118 9.5315 0.0153 0.6605 0.3389 ...
##  $ cd16_cd161                          : num [1:60] 27.1 46.9 55.1 15 23.5 ...
##  $ x25_cd19_cd3_b_cells                : num [1:60] 2.53 8.22 7.86 1.62 9.48 ...
##  $ cd20                                : num [1:60] 0.115 0.111 0.2 0.147 0.102 ...
##  $ cd20_2                              : num [1:60] 2.41 8.09 7.64 1.47 9.38 ...
##  $ x26_cd27_ig_d_naive_b_cells         : num [1:60] 1.19 6.99 4.62 1.29 6.69 ...
##  $ x27_cd27_ig_d_memory_b_cells        : num [1:60] 0.697 0.529 2.22 0.101 0.878 ...
##  $ x28_cd27_ig_d_memory_resting_b_cells: num [1:60] 0.1083 0.2849 0.304 0.0499 1.4241 ...
##  $ x30_cd27_cd38_plasmablasts          : num [1:60] 0.0634 0.0277 0.1582 0.0378 0.062 ...
##  $ cd19_cd3                            : num [1:60] 66.8 35.6 30.2 79.9 65.8 ...
##  $ x31_cd14_monocytes                  : num [1:60] 47.7 19.5 23.2 41.8 48.7 ...
##  $ x32_cd16_non_classical_mono         : num [1:60] 6.24 3.29 3.07 7.15 3.4 ...
##  $ x33_cd16_classical_mono             : num [1:60] 38.5 16.1 19.1 33.5 45.1 ...
##  $ x34_hladr_cd56                      : num [1:60] 43.8 16.7 20.3 33.9 41.1 ...
##  $ cd16_cd123                          : num [1:60] 1.142 0.272 0.657 0.645 0.613 ...
##  $ x35_cd16_cd123_cd11c_p_dc           : num [1:60] 0.709 0.248 0.594 0.37 0.495 ...
##  $ cd16_cd123_2                        : num [1:60] 33.1 10.5 15.1 25.6 37.3 ...
##  $ x36_cd16_cd123_cd11c_m_dc           : num [1:60] 27.7 10.4 14.7 23.6 37.1 ...
##  $ cd123                               : num [1:60] 65.7 35.1 29.5 73.5 65.1 ...
##  $ total_dc                            : num [1:60] 29.5 10.9 15.9 34 38.2 ...
##  $ x37_cd56_cd161_cd123_nk_cells       : num [1:60] 13.01 9.43 2.44 5.33 14.75 ...
##  $ x38_cd16_nk_cells                   : num [1:60] 1.601 1.552 0.874 28.656 5.754 ...
##  $ cd16_nk_cells                       : num [1:60] 11.41 7.85 1.56 41.75 8.9 ...
##  $ cd14_cd11b                          : num [1:60] 46.4 17.3 20.5 35 48.7 ...
##  $ x40_cd14_mdsc_mono                  : num [1:60] 28.983 7.399 9.73 0.215 39.928 ...
##  $ cd66b_cd11b                         : num [1:60] 0.314 0.251 0.265 0.195 0.532 ...
##  $ x41_cd66b_mdsc_grans                : num [1:60] 0.2057 0.0724 0.0985 NA 0.4165 ...
##  $ lutein                              : num [1:60] 735 120 556 419 240 ...
##  $ zeaxanthin                          : num [1:60] 94.6 140.1 99 81.6 53.7 ...
##  $ b_cryptoxanthin                     : num [1:60] 314 286 389 230 570 ...
##  $ a_carotene                          : num [1:60] 195.7 101.8 181.2 112.6 37.1 ...
##  $ b_carotene                          : num [1:60] 927 309 434 536 121 ...
##  $ other_cis_lyc                       : num [1:60] 151 276 161 165 273 ...
##  $ all_trans_lyc                       : num [1:60] 115 259 185 175 336 ...
##  $ x5_cis_lyc                          : num [1:60] 132 262 166 149 235 ...

Wrangle

# convert variables that should be factors to factors
meta_table <- meta_table %>%
  mutate(across(.cols = c("patient_id", "period", 
                          "intervention", "intervention_week", 
                          "pre_post", "sex", "sequence"),
                .fns = as.factor))


# some stuff came in as characters but should be numeric
meta_table <- meta_table %>%
  mutate(across(.cols = c("il_2", "il_10", "il_13", "il_4"),
                .fns = as.numeric))

str(meta_table)

## tibble [60 × 89] (S3: tbl_df/tbl/data.frame)
##  $ patient_id                          : Factor w/ 12 levels "6101","6102",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ sex                                 : Factor w/ 2 levels "F","M": 1 2 1 2 2 2 2 1 2 2 ...
##  $ age_at_enrollment                   : num [1:60] 58 65 60 54 40 57 68 75 36 46 ...
##  $ bmi_at_enrollment                   : num [1:60] 31.1 36.8 30.3 33.1 30.4 ...
##  $ intervention                        : Factor w/ 3 levels "Baseline","Red",..: 2 2 2 3 2 3 2 3 3 3 ...
##  $ sequence                            : Factor w/ 2 levels "R_Y","Y_R": 1 1 1 2 1 2 1 2 2 2 ...
##  $ intervention_week                   : Factor w/ 5 levels "0","2","6","10",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ pre_post                            : Factor w/ 3 levels "baseline","post",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ period                              : Factor w/ 5 levels "B0","B1","B2",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ if_ng                               : num [1:60] 1.16 11.74 0.93 0.63 1.88 ...
##  $ il_1b                               : num [1:60] 15.26 28.75 8.62 11.04 22.5 ...
##  $ il_2                                : num [1:60] 0.98 9.32 0.55 0.65 1.8 0.35 0.38 0.05 NA 0.16 ...
##  $ il_6                                : num [1:60] 2.18 5.83 2.53 1.29 2.25 ...
##  $ il_8                                : num [1:60] 2.86 4.35 3.93 1.93 4.14 3.05 1.35 2.01 4.43 2.86 ...
##  $ il_10                               : num [1:60] 0.05 0.71 1.89 7.17 0.75 2.2 1.91 4.26 0.02 0.48 ...
##  $ il_12p70                            : num [1:60] 2.86 60.93 5.46 2.79 4.06 ...
##  $ mcp_1                               : num [1:60] 51.9 189.3 224.9 244.7 167.9 ...
##  $ tn_fa                               : num [1:60] 14.9 92.3 26.5 19.6 49.6 ...
##  $ il_13                               : num [1:60] 9.64 92.04 21.61 23.17 60.22 ...
##  $ il_5                                : num [1:60] 3 13.04 2.74 5.45 5.38 ...
##  $ il_1ra                              : num [1:60] 3.52 11.36 3.71 3.38 6.08 ...
##  $ il_12p40                            : num [1:60] 51.5 139.4 81 78.7 175.9 ...
##  $ gm_csf                              : num [1:60] 10.9 195.7 17.9 27.4 72.4 ...
##  $ il_4                                : num [1:60] NA 8.8 1.6 4.49 1.22 1.12 0.31 0.91 0.07 NA ...
##  $ x01_cd45_cd66b_lymph_dc_mono        : num [1:60] 99.3 95.6 99.7 99.7 96.4 ...
##  $ x02_cd45_cd66b_grans                : num [1:60] 0.682 0.33 0.305 0.285 0.696 ...
##  $ cd45_cd14                           : num [1:60] 46.9 71 69.8 55.8 42.4 ...
##  $ cd45_cd20                           : num [1:60] 44.2 59.9 60.6 54 31.9 ...
##  $ x03_cd3_cd45_cd3_t_cells            : num [1:60] 25.7 43.7 54.5 14.2 14.3 ...
##  $ x04_tc_rgd_cd3_ab_t_cells           : num [1:60] 22 42.6 53.1 11.2 13.6 ...
##  $ x05_cd4_cd8_cd8_t_cells             : num [1:60] 11.96 3.84 13.23 6.5 5.72 ...
##  $ ccr7_cd8                            : num [1:60] 4.385 1.54 4.186 0.962 1.882 ...
##  $ ccr7_cd8_2                          : num [1:60] 7.56 2.28 9.04 5.52 3.83 ...
##  $ x06_cd45ro_cd45ra_naive_cd8         : num [1:60] 3.692 0.233 3.753 0.652 1.667 ...
##  $ x07_cd46ro_cd45ra_cm_cd8            : num [1:60] 0.506 0.955 0.263 0.196 0.146 ...
##  $ x08_cd45ro_cd45ra_em_cd8            : num [1:60] 6.555 0.395 7.505 4.528 1.807 ...
##  $ x09_cd45r0_cd45ra_te_cd8            : num [1:60] 0.826 1.322 1.2 0.775 1.661 ...
##  $ x10_cd38_hladr_activated_cd8        : num [1:60] 0.2543 0.0823 0.1523 0.1347 0.0742 ...
##  $ x11_cd4_cd8_cd4_t_cells             : num [1:60] 9.14 37.12 37.51 3.64 6.98 ...
##  $ ccr7_cd4                            : num [1:60] 7.98 33.92 32.87 2.49 5.3 ...
##  $ ccr7_cd4_2                          : num [1:60] 1.16 3.19 4.64 1.15 1.68 ...
##  $ x12_cd45ro_cd45ra_naive_cd4         : num [1:60] 3.163 3.77 17.23 0.693 1.731 ...
##  $ x13_cd45ro_cd45ra_cm_cd4            : num [1:60] 3.51 22.39 8.5 1.25 2.63 ...
##  $ x14_cd45ro_cd45ra                   : num [1:60] 1.053 2.964 3.919 0.893 1.488 ...
##  $ x15_cd45ro_cd45ra_te_cd4            : num [1:60] 0.0257 0.0664 0.5457 0.1154 0.092 ...
##  $ x16_cd38_hladr_activated_cd4        : num [1:60] 0.321 0.698 0.384 0.123 0.14 ...
##  $ cd45_cxcr5_th                       : num [1:60] 8.291 28.987 0.133 2.95 6.456 ...
##  $ x17_cd25_cd127_tregs                : num [1:60] 0.348 0.545 0.438 0.131 0.193 ...
##  $ x18_ccr4_cd4_total_ccr4_treg        : num [1:60] 0.317 0.494 0.316 0.123 0.18 ...
##  $ x19_cd45ra_cd45ro_ccr4_treg_naive   : num [1:60] 0.00973 0.03494 0.00676 0.00955 0.00744 ...
##  $ x21_cd45ra_cd45ro_ccr4_treg_memory  : num [1:60] 0.3017 0.4366 0.2905 0.0949 0.166 ...
##  $ x20_hladr_total_ccr4_treg_activated : num [1:60] 0.2867 0.2547 0.1017 0.0838 0.1119 ...
##  $ x22_cxcr3_ccr6_th1                  : num [1:60] 0.7981 2.0038 0.0209 0.4805 0.8362 ...
##  $ x23_cxcr3_ccr6_th2                  : num [1:60] 5.6202 13.4081 0.0747 1.2236 4.3881 ...
##  $ x24_cxcr3_ccr6_th17                 : num [1:60] 1.1118 9.5315 0.0153 0.6605 0.3389 ...
##  $ cd16_cd161                          : num [1:60] 27.1 46.9 55.1 15 23.5 ...
##  $ x25_cd19_cd3_b_cells                : num [1:60] 2.53 8.22 7.86 1.62 9.48 ...
##  $ cd20                                : num [1:60] 0.115 0.111 0.2 0.147 0.102 ...
##  $ cd20_2                              : num [1:60] 2.41 8.09 7.64 1.47 9.38 ...
##  $ x26_cd27_ig_d_naive_b_cells         : num [1:60] 1.19 6.99 4.62 1.29 6.69 ...
##  $ x27_cd27_ig_d_memory_b_cells        : num [1:60] 0.697 0.529 2.22 0.101 0.878 ...
##  $ x28_cd27_ig_d_memory_resting_b_cells: num [1:60] 0.1083 0.2849 0.304 0.0499 1.4241 ...
##  $ x30_cd27_cd38_plasmablasts          : num [1:60] 0.0634 0.0277 0.1582 0.0378 0.062 ...
##  $ cd19_cd3                            : num [1:60] 66.8 35.6 30.2 79.9 65.8 ...
##  $ x31_cd14_monocytes                  : num [1:60] 47.7 19.5 23.2 41.8 48.7 ...
##  $ x32_cd16_non_classical_mono         : num [1:60] 6.24 3.29 3.07 7.15 3.4 ...
##  $ x33_cd16_classical_mono             : num [1:60] 38.5 16.1 19.1 33.5 45.1 ...
##  $ x34_hladr_cd56                      : num [1:60] 43.8 16.7 20.3 33.9 41.1 ...
##  $ cd16_cd123                          : num [1:60] 1.142 0.272 0.657 0.645 0.613 ...
##  $ x35_cd16_cd123_cd11c_p_dc           : num [1:60] 0.709 0.248 0.594 0.37 0.495 ...
##  $ cd16_cd123_2                        : num [1:60] 33.1 10.5 15.1 25.6 37.3 ...
##  $ x36_cd16_cd123_cd11c_m_dc           : num [1:60] 27.7 10.4 14.7 23.6 37.1 ...
##  $ cd123                               : num [1:60] 65.7 35.1 29.5 73.5 65.1 ...
##  $ total_dc                            : num [1:60] 29.5 10.9 15.9 34 38.2 ...
##  $ x37_cd56_cd161_cd123_nk_cells       : num [1:60] 13.01 9.43 2.44 5.33 14.75 ...
##  $ x38_cd16_nk_cells                   : num [1:60] 1.601 1.552 0.874 28.656 5.754 ...
##  $ cd16_nk_cells                       : num [1:60] 11.41 7.85 1.56 41.75 8.9 ...
##  $ cd14_cd11b                          : num [1:60] 46.4 17.3 20.5 35 48.7 ...
##  $ x40_cd14_mdsc_mono                  : num [1:60] 28.983 7.399 9.73 0.215 39.928 ...
##  $ cd66b_cd11b                         : num [1:60] 0.314 0.251 0.265 0.195 0.532 ...
##  $ x41_cd66b_mdsc_grans                : num [1:60] 0.2057 0.0724 0.0985 NA 0.4165 ...
##  $ lutein                              : num [1:60] 735 120 556 419 240 ...
##  $ zeaxanthin                          : num [1:60] 94.6 140.1 99 81.6 53.7 ...
##  $ b_cryptoxanthin                     : num [1:60] 314 286 389 230 570 ...
##  $ a_carotene                          : num [1:60] 195.7 101.8 181.2 112.6 37.1 ...
##  $ b_carotene                          : num [1:60] 927 309 434 536 121 ...
##  $ other_cis_lyc                       : num [1:60] 151 276 161 165 273 ...
##  $ all_trans_lyc                       : num [1:60] 115 259 185 175 336 ...
##  $ x5_cis_lyc                          : num [1:60] 132 262 166 149 235 ...

# changing factor levels for pre_post
meta_table$pre_post <- factor(meta_table$pre_post,
                              levels = c("pre", "post"))

levels(meta_table$pre_post)  

## [1] "pre"  "post"

# Calculate total_cis_lyc, total_lyc, and total_carotenoids
meta_table <- meta_table %>%
  rename(n5_cis_lyc = x5_cis_lyc) %>%
  mutate(total_cis_lyc = other_cis_lyc + n5_cis_lyc,
         total_lyc = all_trans_lyc + total_cis_lyc,
         total_carotenoids = lutein + zeaxanthin + b_cryptoxanthin + 
                             a_carotene + b_carotene + total_lyc) 

Carotenoids

All-trans-lyc levels

# line plots for each subject at each timepoint
meta_table %>% 
  ggplot(aes(x = intervention_week, y = all_trans_lyc, color = intervention)) +
  geom_point() + 
  geom_line(aes(group = intervention)) +
  scale_color_manual(values = c("Baseline" = "gray", 
                                           "Yellow" = "gold",
                                           "Red" = "tomato1")) +
  facet_wrap(vars(patient_id), scales = "free_y") + 
  theme_bw() +
  labs(x = "Intervention Week",
       y = "All-trans-lycopene levels (nmol/L)",
       title = "All-trans-lycopene levels in each patient before/after each intervention")

Total cis-lyc levels

# line plots for each subject at each timepoint
meta_table %>% 
  ggplot(aes(x = intervention_week, y = total_cis_lyc, color = intervention)) +
  geom_point() + 
  geom_line(aes(group = intervention)) +
  scale_color_manual(values = c("Baseline" = "gray", 
                                           "Yellow" = "gold",
                                           "Red" = "tomato1")) +
  facet_wrap(vars(patient_id)) +
  theme_bw() +
  labs(x = "Intervention Week",
       y = "Total cis-lycopene levels (nmol/L)",
       title = "Total cis-lycopene levels in each patient before/after each intervention")

Total lyc levels

lineplots

# line plots for each subject at each timepoint
meta_table %>% 
  ggplot(aes(x = intervention_week, y = total_lyc, color = intervention)) +
  geom_point() + 
  geom_line(aes(group = intervention)) +
  scale_color_manual(values = c("Baseline" = "gray", 
                                           "Yellow" = "gold",
                                           "Red" = "tomato1")) +
  facet_wrap(vars(patient_id), scales = "free_y") +
  theme_bw() +
  labs(x = "Intervention Week",
       y = "Total lycopene levels (nmol/L)",
       title = "Total lycopene levels in each patient before/after each intervention")

Boxplots

wrangling

# create a more specific pre_post_intervention column
meta_table_edited <- meta_table %>%
  unite(col = "pre_post_intervention",
        c("pre_post","intervention"),
        sep = "_",
        remove = FALSE)

# make pre_post_intervention column factors
meta_table_edited$pre_post_intervention <- as.factor(meta_table_edited$pre_post_intervention)

# relevel factor columns
meta_table_edited$pre_post_intervention <- factor(meta_table_edited$pre_post_intervention, levels = c("pre_Yellow", "post_Yellow", "pre_Red", "post_Red"))

meta_table_edited$intervention <- factor(meta_table_edited$intervention,
                                         levels = c("Yellow", "Red"))

# make legend title
legendtitle_ppintervention <- "Timepoint"


# labels
labs_ppintervention <- c("pre control",
                         "post control",
                         "pre Tomato-Soy",
                         "post Tomato-Soy")

figure

meta_table_edited %>% 
  filter(intervention != "Baseline") %>%
  ggplot(aes(x = intervention, y = total_lyc, fill = pre_post_intervention)) +
  geom_boxplot(outlier.shape = NA) + 
  scale_fill_manual(legendtitle_ppintervention,
                    values = c("pre_Red" = "#FF9966",
                               "post_Red" = "#FF3300",
                               "pre_Yellow" = "#FFFF99",
                               "post_Yellow" = "yellow"),
                    labels = labs_ppintervention) +
  theme_clean() +
  labs(x = "",
       y = "Total lycopene levels (nmol/L)",
       title = "Total lycopene levels before/after juice interventions")

meta_table_edited %>% 
  filter(intervention != "Baseline") %>%
  ggplot(aes(x = pre_post_intervention, y = total_lyc, fill = intervention)) +
  geom_boxplot(outlier.shape = NA, show.legend = FALSE) + 
  scale_fill_manual(values = c("Red" = "tomato1",
                               "Yellow" = "yellow1")) +
  scale_x_discrete(labels = labs_ppintervention) +
  theme_clean(base_size = 22, base_family = "sans") +
  labs(x = "",
       y = "Overall plasma lycopene conc. (nmol/L)",
       title = "Total plasma lycopene levels",
       subtitle = "Before and after juice interventions")

# publish-ready plot
library(ggpubr)

(total_lyc_levels <- meta_table_edited %>% 
    filter(intervention != "Baseline") %>%
  ggpaired(x = "pre_post", y = "total_lyc", fill = "intervention", line.color = "gray", line.size = 1, facet.by = "intervention", short.panel.labs = FALSE, panel.labs = list(intervention = c("", ""))) +
  scale_fill_manual(values = c("Red" = "tomato1",
                               "Yellow" = "yellow1"),
                    labels = c("Control", "Tomato-Soy"),
                    name = "Intervention") +
  geom_line(aes(group = patient_id), colour = "gray", linewidth = 0.15) +
  theme_clean(base_size = 18, base_family = "sans") +
  labs(x = "",
       y = "nmol/L plasma",
       title = "Concentration of Lycopene",
       subtitle = ""))

export

ggsave(filename = "/Users/mariasholola/Documents/GitHub/USDA-Inflammation-Metabolomics/Plots/total-lyc-red-and-yellow-boxplots.svg", plot = total_lyc_levels, width = 12)

Descriptive stats

wrangle

# convert meta_table_edited to long format for lycopene
meta_table_lyc_long <- meta_table_edited %>%
  pivot_longer(cols = total_lyc,
               names_to = "total_lycopene",
               values_to = "nmol_per_L")

Avg/std dev

tomsoy

meta_table_lyc_long %>%
  filter(intervention == "Red") %>%
  group_by(pre_post) %>%
  summarize(mean = mean(nmol_per_L),
            stdev = sd(nmol_per_L))

## # A tibble: 2 × 3
##   pre_post  mean stdev
##   <fct>    <dbl> <dbl>
## 1 pre       523.  228.
## 2 post     1298.  665.

yellow

meta_table_lyc_long %>%
  filter(intervention == "Yellow") %>%
  group_by(pre_post) %>%
  summarize(mean = mean(nmol_per_L),
            stdev = sd(nmol_per_L))

## # A tibble: 2 × 3
##   pre_post  mean stdev
##   <fct>    <dbl> <dbl>
## 1 pre       721.  377.
## 2 post      703.  434.

Mean fold changes

lyc_long_subset <- meta_table_lyc_long %>%
  select(patient_id, nmol_per_L, pre_post_intervention) %>%
  pivot_wider(names_from = pre_post_intervention,
              values_from = nmol_per_L) %>%
  mutate(red_FC = post_Red/pre_Red,
         yellow_FC = post_Yellow/pre_Yellow)

lyc_long_subset %>%
  summarize(mean_red_FC = mean(red_FC),
            mean_yellow_FC = mean(yellow_FC))

## # A tibble: 1 × 2
##   mean_red_FC mean_yellow_FC
##         <dbl>          <dbl>
## 1        2.83          0.966

Normality checks

Plot histogram

gghistogram(meta_table_lyc_long$nmol_per_L, bins = 40)

Shapiro’s normality test

# shapiro normality test for total lycopene 

meta_table_lyc_long %>%
  group_by(intervention) %>%
  shapiro_test(vars = "nmol_per_L")

## # A tibble: 3 × 4
##   intervention variable   statistic       p
##   <fct>        <chr>          <dbl>   <dbl>
## 1 Yellow       nmol_per_L     0.876 0.00698
## 2 Red          nmol_per_L     0.848 0.00202
## 3 <NA>         nmol_per_L     0.913 0.233

P val for shapiro test turned out to < 0.05 for both interventions, suggesting data here is not normal. Let’s log transform before running paired t-tests.

Log transform

meta_table_lyc_log2_long <- meta_table_lyc_long %>%
  mutate(log2_nmol_per_L = log2(nmol_per_L))

Without transforming, here’s a paired t-test

compare_means(nmol_per_L ~ pre_post, meta_table_lyc_log2_long, method = "t.test", paired = TRUE, group.by = "intervention")

## # A tibble: 2 × 9
##   intervention .y.        group1 group2       p  p.adj p.format p.signif method
##   <fct>        <chr>      <chr>  <chr>    <dbl>  <dbl> <chr>    <chr>    <chr> 
## 1 Red          nmol_per_L pre    post   0.00104 0.0021 0.001    **       T-test
## 2 Yellow       nmol_per_L pre    post   0.709   0.71   0.709    ns       T-test

normality check

Plot histogram

gghistogram(meta_table_lyc_log2_long$log2_nmol_per_L, bins = 40)

Shapiro’s normality test

# shapiro normality test for total lycopene 
meta_table_lyc_log2_long %>%
  group_by(intervention) %>%
  shapiro_test(vars = "log2_nmol_per_L")

## # A tibble: 3 × 4
##   intervention variable        statistic      p
##   <fct>        <chr>               <dbl>  <dbl>
## 1 Yellow       log2_nmol_per_L     0.928 0.0869
## 2 Red          log2_nmol_per_L     0.964 0.521 
## 3 <NA>         log2_nmol_per_L     0.903 0.173

paired t test

compare_means(log2_nmol_per_L ~ pre_post, meta_table_lyc_log2_long, method = "t.test", paired = TRUE, group.by = "intervention")

## # A tibble: 2 × 9
##   intervention .y.        group1 group2       p   p.adj p.format p.signif method
##   <fct>        <chr>      <chr>  <chr>    <dbl>   <dbl> <chr>    <chr>    <chr> 
## 1 Red          log2_nmol… pre    post   3.84e-5 7.70e-5 3.8e-05  ****     T-test
## 2 Yellow       log2_nmol… pre    post   3.21e-1 3.2 e-1 0.32     ns       T-test

Total lycopene levels are significantly increasing only after post-Red intervention

Cytokines

# convert cytokines from wide to long
cytokines_long <- meta_table[-c(25:81)] %>%
  pivot_longer(cols = if_ng:il_4,
               names_to = "cytokines",
               values_to = "cyto_conc_pg_ml")

#cytokines_long_nona <- na.omit(cytokines_long)

note 11/30/23: I took this part out since NAs don’t matter for lmm

Mixed linear modeling

After testing several models, the best model turned out to have pre_post as a fixed variable and patient_id as random effect. Carotenoids added no statistically significant effect to the model, suggesting that carotenoids are not contributing to cytokine levels. Set REML = FALSE when comparing models.

Sequence effect test

# turn off sci notation
options(scipen = 999)

# sequence effect test lmer model function
cyto_seq_function <- function(cytokines_long) {
  lmer(cyto_conc_pg_ml ~ sequence + (1|patient_id), 
       data = cytokines_long, 
       REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_seq_models <- map(split(cytokines_long,
                         cytokines_long$cytokines),
                   cyto_seq_function)

# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_seq_results <- map_df(cyto_seq_models,
                          tidy,
                          .id = "cytokines")

# extract fixed coefficients for sequence only
cyto_seq_results_coefonly <- cyto_seq_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

# results 
kable(cyto_seq_results_coefonly, format = "markdown", digits = 3)

cytokines	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
gm_csf	fixed	NA	sequenceY_R	-27.444	24.681	-1.112	10.000	0.292	ns
if_ng	fixed	NA	sequenceY_R	-1.194	1.477	-0.808	10.000	0.438	ns
il_10	fixed	NA	sequenceY_R	1.728	1.220	1.416	9.996	0.187	ns
il_12p40	fixed	NA	sequenceY_R	-16.115	23.825	-0.676	10.000	0.514	ns
il_12p70	fixed	NA	sequenceY_R	-7.289	7.609	-0.958	10.000	0.361	ns
il_13	fixed	NA	sequenceY_R	-18.773	13.629	-1.377	10.724	0.196	ns
il_1b	fixed	NA	sequenceY_R	-1.552	8.358	-0.186	10.000	0.856	ns
il_1ra	fixed	NA	sequenceY_R	-2.149	1.685	-1.275	10.000	0.231	ns
il_2	fixed	NA	sequenceY_R	-1.585	1.440	-1.101	8.191	0.302	ns
il_4	fixed	NA	sequenceY_R	-0.667	1.778	-0.375	6.000	0.721	ns
il_5	fixed	NA	sequenceY_R	1.657	3.027	0.547	10.000	0.596	ns
il_6	fixed	NA	sequenceY_R	-2.432	1.896	-1.283	10.000	0.229	ns
il_8	fixed	NA	sequenceY_R	0.432	0.674	0.641	10.000	0.536	ns
mcp_1	fixed	NA	sequenceY_R	49.779	28.416	1.752	10.000	0.110	ns
tn_fa	fixed	NA	sequenceY_R	-6.730	10.112	-0.666	10.000	0.521	ns

Are there any significant cytokines?

# extract statistically significant cytokines 
cytokines_seq_psig <- cyto_seq_results_coefonly %>%
  filter(cyto_seq_results_coefonly$p.value < 0.05)

print(cytokines_seq_psig$cytokines)

## character(0)

• No statistically significant cytokines, suggesting there are no sequence effects

Intervention comparison

# filter data for post-intervention only
cytokines_post_long <- cytokines_long %>%
  filter(intervention != "Baseline") %>%
  filter(pre_post == "post")

# treatment effect test lmer model function
cyto_trt_function <- function(cytokines_post_long) {
  lmer(cyto_conc_pg_ml ~ intervention + (1|patient_id), 
       data = cytokines_post_long, 
       REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_trt_models <- map(split(cytokines_post_long,
                         cytokines_post_long$cytokines),
                   cyto_trt_function)

# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_trt_results <- map_df(cyto_trt_models,
                          tidy,
                          .id = "cytokines")

# extract fixed coefficients for sequence only
cyto_trt_coefonly <- cyto_trt_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Results

kable(cyto_trt_coefonly, format = "markdown", digits = 3)

cytokines	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
gm_csf	fixed	NA	interventionYellow	-2.237	2.904	-0.770	11.000	0.457	ns
if_ng	fixed	NA	interventionYellow	0.081	0.213	0.380	11.000	0.711	ns
il_10	fixed	NA	interventionYellow	0.121	0.158	0.766	10.035	0.461	ns
il_12p40	fixed	NA	interventionYellow	6.302	5.595	1.126	11.000	0.284	ns
il_12p70	fixed	NA	interventionYellow	-0.220	0.361	-0.609	11.000	0.555	ns
il_13	fixed	NA	interventionYellow	-0.213	1.716	-0.124	7.124	0.905	ns
il_1b	fixed	NA	interventionYellow	1.167	1.119	1.043	11.000	0.319	ns
il_1ra	fixed	NA	interventionYellow	0.311	0.469	0.662	11.000	0.521	ns
il_2	fixed	NA	interventionYellow	-0.153	0.168	-0.909	6.131	0.398	ns
il_4	fixed	NA	interventionYellow	-0.129	0.075	-1.713	7.000	0.131	ns
il_5	fixed	NA	interventionYellow	0.500	0.708	0.706	11.000	0.495	ns
il_6	fixed	NA	interventionYellow	0.991	0.750	1.322	11.000	0.213	ns
il_8	fixed	NA	interventionYellow	-0.060	0.355	-0.169	11.000	0.869	ns
mcp_1	fixed	NA	interventionYellow	3.200	8.404	0.381	11.000	0.711	ns
tn_fa	fixed	NA	interventionYellow	0.458	1.097	0.418	11.000	0.684	ns

# extract statistically significant cytokines 
cytokines_trt_psig <- cyto_trt_coefonly %>%
  filter(cyto_trt_coefonly$p.value < 0.05)


print(cytokines_trt_psig$cytokines)

## character(0)

• There are no cytokines with significantly different levels when comparing post-red vs. post-yellow treatment

Yellow treatment effects

# subset data into yellow results
cytokines_Y <- cytokines_long %>%
  filter(intervention == "Yellow")

# treatment effect test lmer model function
cyto_trtY_function <- function(cytokines_Y) {
  lmer(cyto_conc_pg_ml ~ pre_post + (1|patient_id), 
       data = cytokines_Y, 
       REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_trtY_models <- map(split(cytokines_Y,
                         cytokines_Y$cytokines),
                   cyto_trtY_function)

# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_trtY_results <- map_df(cyto_trtY_models,
                          tidy,
                          .id = "cytokines")

# extract fixed coefficients for sequence only
cyto_trtY_coefonly <- cyto_trtY_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Results

kable(cyto_trtY_coefonly, format = "markdown", digits = 3)

cytokines	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
gm_csf	fixed	NA	pre_postpost	-10.773	9.940	-1.084	11.000	0.302	ns
if_ng	fixed	NA	pre_postpost	-0.582	0.554	-1.050	11.000	0.316	ns
il_10	fixed	NA	pre_postpost	0.003	0.144	0.024	9.060	0.982	ns
il_12p40	fixed	NA	pre_postpost	-7.124	9.540	-0.747	11.000	0.471	ns
il_12p70	fixed	NA	pre_postpost	-3.552	3.510	-1.012	11.000	0.333	ns
il_13	fixed	NA	pre_postpost	-10.171	8.663	-1.174	7.959	0.274	ns
il_1b	fixed	NA	pre_postpost	-4.913	3.071	-1.600	11.000	0.138	ns
il_1ra	fixed	NA	pre_postpost	-1.743	1.465	-1.189	11.000	0.259	ns
il_2	fixed	NA	pre_postpost	-1.097	0.952	-1.151	7.491	0.285	ns
il_4	fixed	NA	pre_postpost	-0.752	0.722	-1.043	7.000	0.332	ns
il_5	fixed	NA	pre_postpost	-0.906	1.413	-0.641	11.000	0.535	ns
il_6	fixed	NA	pre_postpost	-0.614	0.659	-0.932	11.000	0.372	ns
il_8	fixed	NA	pre_postpost	0.130	0.296	0.439	11.000	0.669	ns
mcp_1	fixed	NA	pre_postpost	1.432	5.523	0.259	11.000	0.800	ns
tn_fa	fixed	NA	pre_postpost	-3.615	4.206	-0.859	11.000	0.408	ns

Significant cytokines

# extract statistically significant cytokines 
cytokines_trtY_psig <- cyto_trtY_coefonly %>%
  filter(cyto_trtY_coefonly$p.value < 0.05)


print(cytokines_trtY_psig$cytokines)

## character(0)

• No significant cytokines within yellow treatment

Red treatment effects

# subset data into yellow results
cytokines_R <- cytokines_long %>%
  filter(intervention == "Red")

# treatment effect test lmer model function
cyto_trtR_function <- function(cytokines_R) {
  lmer(cyto_conc_pg_ml ~ pre_post + (1|patient_id), 
       data = cytokines_R, 
       REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_trtR_models <- map(split(cytokines_R,
                         cytokines_R$cytokines),
                   cyto_trtR_function)

# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_trtR_results <- map_df(cyto_trtR_models,
                          tidy,
                          .id = "cytokines")

# extract fixed coefficients for sequence only
cyto_trtR_coefonly <- cyto_trtR_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Results

kable(cyto_trtR_coefonly, format = "markdown", digits = 3) %>%
  row_spec(11, bold = TRUE)

cytokines	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
gm_csf	fixed	NA	pre_postpost	-7.572	5.759	-1.315	11.000	0.215	ns
if_ng	fixed	NA	pre_postpost	-0.497	0.412	-1.209	11.000	0.252	ns
il_10	fixed	NA	pre_postpost	-0.157	0.139	-1.131	10.028	0.284	ns
il_12p40	fixed	NA	pre_postpost	-13.373	6.652	-2.010	11.000	0.070	ns
il_12p70	fixed	NA	pre_postpost	-3.119	2.415	-1.291	11.000	0.223	ns
il_13	fixed	NA	pre_postpost	-5.090	4.748	-1.072	7.669	0.316	ns
il_1b	fixed	NA	pre_postpost	-1.908	1.306	-1.461	11.000	0.172	ns
il_1ra	fixed	NA	pre_postpost	-0.435	0.453	-0.961	11.000	0.357	ns
il_2	fixed	NA	pre_postpost	-0.522	0.556	-0.938	7.313	0.378	ns
il_4	fixed	NA	pre_postpost	-0.581	0.527	-1.104	7.000	0.306	ns
il_5	fixed	NA	pre_postpost	-1.808	0.815	-2.220	11.000	0.048	*
il_6	fixed	NA	pre_postpost	-1.102	0.820	-1.343	11.000	0.206	ns
il_8	fixed	NA	pre_postpost	0.073	0.217	0.338	11.000	0.741	ns
mcp_1	fixed	NA	pre_postpost	-2.799	6.896	-0.406	11.000	0.693	ns
tn_fa	fixed	NA	pre_postpost	-4.860	2.934	-1.656	11.000	0.126	ns

Significant cytokines

# extract statistically significant cytokines 
cytokines_trtR_psig <- cyto_trtR_coefonly %>%
  filter(cyto_trtR_coefonly$p.value < 0.05)

print(cytokines_trtR_psig$cytokines)

## [1] "il_5"

• IL-5 is the only significantly different cytokine within red intervention. Same p value from t-test (= 0.0484)
• IL-5 on average goes down by ~1.8 pg/mL from pre to post in red intervention

IL-5 plots

lineplots

meta_table %>%
  filter(intervention == "Red") %>%
  ggplot(aes(x = intervention_week, y = il_5, 
             color = intervention, group = intervention)) +
  geom_line() +
  theme_classic() +
  labs(x = "Week",
       y = "IL-5 levels (pg/mL)",
       title = "IL-5 changes by subject before and after red intervention",
       color = "Patient ID") +
  facet_wrap(~patient_id,
             scales = "free_y") +
  scale_color_manual(values = c("Red" = "tomato1"))

meta_table %>% 
  filter(intervention == "Red") %>%
  ggplot(aes(x = pre_post, y = il_5, color = patient_id)) +
  geom_line(aes(group = patient_id)) +
  labs(x = "",
       y = "IL-5 conc (pg/mL)",
       title = "IL-5 levels in each patient before and after red intervention") +
  theme_classic()

meta_table %>% 
  filter(intervention != "Baseline") %>%
  ggplot(aes(x = pre_post, y = il_5, color = intervention)) +
  geom_line(aes(group = intervention)) +
  scale_color_manual(values = c("Yellow" = "gold",
                                "Red" = "tomato1")) +
  facet_wrap(vars(patient_id), scales = "free_y") + 
  theme_classic() +
  labs(x = "",
       y = "IL-5 conc (pg/mL)",
       title = "IL-5 levels in each patient pre- and post- red and yellow intervention")

boxplots

cytokines_long %>%
  filter(intervention == "Red") %>%
  filter(cytokines == "il_5") %>%
  ggplot(aes(x = pre_post, y = cyto_conc_pg_ml, fill = intervention)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter() +
  scale_fill_manual(values = c("Red" = "tomato1")) +
  labs(x = "",
       y = "IL- 5 levels (pg/mL)",
       title = "IL-5 levels before and after red intervention") +
  theme_bw()

• IL-5 on average goes down by ~1.8 pg/mL from pre to post in red intervention. This is a significant decrease according to the mixed linear model.

Log-transformed cytokines

I know the data isn’t normally distributed. Mixed linear models do not assume normality of the dependant variable, however the residuals should be normally distributed.

resid normality checks

I need to figure out how to make a function that extracts makes a QQplot for each model. But since IL-5 was significant, let’s look at how the residuals look for this cytokine (for every comparison) before log transforming

#qqplot_fx <- function(cytokine){
  #qqnorm(resid(cyto_seq_models$cytokine))
  #qqline(resid(cyto_seq_models$cytokine))
#}

#qqplot_fx("il_5")

#cyto_seq_results$

seq

qqnorm(resid(cyto_seq_models$il_5))
qqline(resid(cyto_seq_models$il_5))

That doesn’t look good.

intervention

qqnorm(resid(cyto_trt_models$il_5))
qqline(resid(cyto_trt_models$il_5))

doesn’t look as bad as sequence

red

qqnorm(resid(cyto_trtR_models$il_5))
qqline(resid(cyto_trtR_models$il_5))

doesn’t look that bad for Red

yellow

qqnorm(resid(cyto_trtY_models$il_5))
qqline(resid(cyto_trtY_models$il_5))

I will log transform the data and perform mixed linear modeling and paired t-tests for every comparison to see how the data compares.

Wrangle

cytokines_log2conc_long <- cytokines_long %>%
  mutate(log2_cyto_conc_pg_ml = log2(cyto_conc_pg_ml))

Mixed linear modeling

Sequence effect test

We want to make sure there are no seq effects

# sequence effect test lmer model function
cyto_log2_seq_function <- function(cytokines_log2conc_long) {
  lmer(log2_cyto_conc_pg_ml ~ sequence + (1|patient_id), 
       data = cytokines_log2conc_long, 
       REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_log2_seq_models <- map(split(cytokines_log2conc_long,
                         cytokines_log2conc_long$cytokines),
                   cyto_log2_seq_function)

# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_log2_seq_results <- map_df(cyto_log2_seq_models,
                          tidy,
                          .id = "cytokines")

# extract fixed coefficients for sequence only
cyto_log2_seq_results_coefonly <- cyto_log2_seq_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Resid check

# before log transform
qqnorm(resid(cyto_seq_models$il_5))
qqline(resid(cyto_seq_models$il_5))

# after log transform
qqnorm(resid(cyto_log2_seq_models$il_5))
qqline(resid(cyto_log2_seq_models$il_5))

This looks better after log transform.

Results

# results 
kable(cyto_log2_seq_results_coefonly, format = "markdown", digits = 3)

cytokines	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
gm_csf	fixed	NA	sequenceY_R	-0.703	0.813	-0.864	10.000	0.408	ns
if_ng	fixed	NA	sequenceY_R	-0.639	1.065	-0.600	10.000	0.562	ns
il_10	fixed	NA	sequenceY_R	0.897	1.380	0.650	9.855	0.530	ns
il_12p40	fixed	NA	sequenceY_R	-0.181	0.408	-0.444	10.000	0.666	ns
il_12p70	fixed	NA	sequenceY_R	-0.804	0.889	-0.905	10.000	0.387	ns
il_13	fixed	NA	sequenceY_R	-1.081	0.848	-1.275	10.230	0.230	ns
il_1b	fixed	NA	sequenceY_R	-0.321	0.836	-0.384	10.000	0.709	ns
il_1ra	fixed	NA	sequenceY_R	-0.570	0.518	-1.101	10.000	0.297	ns
il_2	fixed	NA	sequenceY_R	-1.546	1.196	-1.292	8.035	0.232	ns
il_4	fixed	NA	sequenceY_R	-0.557	1.491	-0.374	6.000	0.721	ns
il_5	fixed	NA	sequenceY_R	0.115	0.595	0.194	10.000	0.850	ns
il_6	fixed	NA	sequenceY_R	-0.976	0.607	-1.607	10.000	0.139	ns
il_8	fixed	NA	sequenceY_R	0.287	0.371	0.773	10.000	0.457	ns
mcp_1	fixed	NA	sequenceY_R	0.469	0.304	1.543	10.000	0.154	ns
tn_fa	fixed	NA	sequenceY_R	-0.033	0.404	-0.083	10.000	0.936	ns

Are there any significant cytokines?

# extract statistically significant cytokines 
cytokines_log2_seq_psig <- cyto_log2_seq_results_coefonly %>%
  filter(cyto_log2_seq_results_coefonly$p.value < 0.05)

print(cytokines_log2_seq_psig$cytokines)

## character(0)

• No statistically significant cytokines, suggesting there are no sequence effects

Intervention comparison

# filter data for post-intervention only
cytokines_log2conc_post_long <- cytokines_log2conc_long %>%
  filter(intervention != "Baseline") %>%
  filter(pre_post == "post")

# treatment effect test lmer model function
cyto_log2_trt_function <- function(cytokines_log2conc_post_long) {
  lmer(log2_cyto_conc_pg_ml ~ intervention + (1|patient_id), 
       data = cytokines_log2conc_post_long, 
       REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_log2_trt_models <- map(split(cytokines_log2conc_post_long,
                         cytokines_log2conc_post_long$cytokines),
                   cyto_log2_trt_function)

# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_log2_trt_results <- map_df(cyto_log2_trt_models,
                          tidy,
                          .id = "cytokines")

# extract fixed coefficients for sequence only
cyto_log2_trt_coefonly <- cyto_log2_trt_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Resid check

# before log transform
qqnorm(resid(cyto_trt_models$il_5))
qqline(resid(cyto_trt_models$il_5))

# log transformed
qqnorm(resid(cyto_log2_trt_models$il_5))
qqline(resid(cyto_log2_trt_models$il_5))

Don’t know if this helps for intervention comparison

Results

kable(cyto_log2_trt_coefonly, format = "markdown", digits = 3)

cytokines	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
gm_csf	fixed	NA	interventionYellow	0.116	0.118	0.980	11.000	0.348	ns
if_ng	fixed	NA	interventionYellow	0.161	0.230	0.700	11.000	0.498	ns
il_10	fixed	NA	interventionYellow	0.183	0.110	1.665	9.985	0.127	ns
il_12p40	fixed	NA	interventionYellow	0.177	0.153	1.159	11.000	0.271	ns
il_12p70	fixed	NA	interventionYellow	0.166	0.153	1.083	11.000	0.302	ns
il_13	fixed	NA	interventionYellow	0.146	0.146	0.997	7.041	0.352	ns
il_1b	fixed	NA	interventionYellow	0.126	0.164	0.767	11.000	0.459	ns
il_1ra	fixed	NA	interventionYellow	0.083	0.155	0.535	11.000	0.603	ns
il_2	fixed	NA	interventionYellow	0.216	0.291	0.745	6.079	0.484	ns
il_4	fixed	NA	interventionYellow	-0.052	0.070	-0.743	7.000	0.482	ns
il_5	fixed	NA	interventionYellow	0.379	0.216	1.749	11.000	0.108	ns
il_6	fixed	NA	interventionYellow	0.331	0.263	1.260	11.000	0.234	ns
il_8	fixed	NA	interventionYellow	-0.105	0.156	-0.674	11.000	0.514	ns
mcp_1	fixed	NA	interventionYellow	0.007	0.066	0.108	11.000	0.916	ns
tn_fa	fixed	NA	interventionYellow	0.035	0.071	0.499	11.000	0.627	ns

# extract statistically significant cytokines 
cytokines_log2_trt_psig <- cyto_log2_trt_coefonly %>%
  filter(cyto_log2_trt_coefonly$p.value < 0.05)


print(cytokines_log2_trt_psig$cytokines)

## character(0)

• There are no cytokines with significantly different levels when comparing post-red vs. post-yellow treatment

Yellow treatment effects

# subset data into yellow results
cytokines_Y_log2 <- cytokines_log2conc_long %>%
  filter(intervention == "Yellow")

# treatment effect test lmer model function
cyto_trtY_log2_function <- function(cytokines_Y_log2) {
  lmer(log2_cyto_conc_pg_ml ~ pre_post + (1|patient_id), 
       data = cytokines_Y_log2, 
       REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_trtY_log2_models <- map(split(cytokines_Y_log2,
                         cytokines_Y_log2$cytokines),
                   cyto_trtY_log2_function)

# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_trtY_log2_results <- map_df(cyto_trtY_log2_models,
                          tidy,
                          .id = "cytokines")

# extract fixed coefficients for sequence only
cyto_trtY_log2_coefonly <- cyto_trtY_log2_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Resid check

# before log transform 
qqnorm(resid(cyto_trtY_models$il_5))
qqline(resid(cyto_trtY_models$il_5))

# after log transform
qqnorm(resid(cyto_trtY_log2_models$il_5))
qqline(resid(cyto_trtY_log2_models$il_5))

Results

kable(cyto_trtY_log2_coefonly, format = "markdown", digits = 3)

cytokines	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
gm_csf	fixed	NA	pre_postpost	-0.183	0.135	-1.353	11.000	0.203	ns
if_ng	fixed	NA	pre_postpost	-0.235	0.165	-1.423	11.000	0.182	ns
il_10	fixed	NA	pre_postpost	-0.241	0.226	-1.066	8.920	0.314	ns
il_12p40	fixed	NA	pre_postpost	-0.068	0.126	-0.543	11.000	0.598	ns
il_12p70	fixed	NA	pre_postpost	-0.053	0.145	-0.365	11.000	0.722	ns
il_13	fixed	NA	pre_postpost	-0.130	0.242	-0.539	6.933	0.607	ns
il_1b	fixed	NA	pre_postpost	-0.206	0.208	-0.991	11.000	0.343	ns
il_1ra	fixed	NA	pre_postpost	-0.213	0.220	-0.970	11.000	0.353	ns
il_2	fixed	NA	pre_postpost	-0.029	0.338	-0.087	7.055	0.933	ns
il_4	fixed	NA	pre_postpost	-0.121	0.192	-0.633	7.000	0.547	ns
il_5	fixed	NA	pre_postpost	0.124	0.239	0.518	11.000	0.615	ns
il_6	fixed	NA	pre_postpost	-0.125	0.356	-0.350	11.000	0.733	ns
il_8	fixed	NA	pre_postpost	-0.024	0.141	-0.170	11.000	0.868	ns
mcp_1	fixed	NA	pre_postpost	-0.012	0.043	-0.278	11.000	0.786	ns
tn_fa	fixed	NA	pre_postpost	-0.082	0.094	-0.876	11.000	0.400	ns

Significant cytokines

# extract statistically significant cytokines 
cytokines_trtY_log2_psig <- cyto_trtY_log2_coefonly %>%
  filter(cyto_trtY_log2_coefonly$p.value < 0.05)


print(cytokines_trtY_log2_psig$cytokines)

## character(0)

• No significant cytokines within yellow treatment

Red treatment effects

# subset data into yellow results
cytokines_R_log2 <- cytokines_log2conc_long %>%
  filter(intervention == "Red")

# treatment effect test lmer model function
cyto_trtR_log2_function <- function(cytokines_R_log2) {
  lmer(log2_cyto_conc_pg_ml ~ pre_post + (1|patient_id), 
       data = cytokines_R_log2, 
       REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_trtR_log2_models <- map(split(cytokines_R_log2,
                         cytokines_R_log2$cytokines),
                   cyto_trtR_log2_function)

# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_trtR_log2_results <- map_df(cyto_trtR_log2_models,
                          tidy,
                          .id = "cytokines")

# extract fixed coefficients for sequence only
cyto_trtR_log2coefonly <- cyto_trtR_log2_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Resid check

# before log transform
qqnorm(resid(cyto_trtR_models$il_5))
qqline(resid(cyto_trtR_models$il_5))

# log transformed
qqnorm(resid(cyto_trtR_log2_models$il_5))
qqline(resid(cyto_trtR_log2_models$il_5))

Results

kable(cyto_trtR_log2coefonly, format = "markdown", digits = 3) 

cytokines	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
gm_csf	fixed	NA	pre_postpost	-0.308	0.103	-2.985	11.000	0.012	*
if_ng	fixed	NA	pre_postpost	-0.307	0.197	-1.559	11.000	0.147	ns
il_10	fixed	NA	pre_postpost	-0.137	0.422	-0.325	9.926	0.752	ns
il_12p40	fixed	NA	pre_postpost	-0.288	0.156	-1.840	11.000	0.093	ns
il_12p70	fixed	NA	pre_postpost	-0.622	0.240	-2.596	11.000	0.025	*
il_13	fixed	NA	pre_postpost	-0.130	0.199	-0.652	7.181	0.535	ns
il_1b	fixed	NA	pre_postpost	-0.227	0.164	-1.387	11.000	0.193	ns
il_1ra	fixed	NA	pre_postpost	-0.109	0.120	-0.912	11.000	0.381	ns
il_2	fixed	NA	pre_postpost	-0.273	0.277	-0.988	7.142	0.356	ns
il_4	fixed	NA	pre_postpost	-0.200	0.156	-1.282	7.000	0.241	ns
il_5	fixed	NA	pre_postpost	-0.373	0.172	-2.166	11.000	0.053	ns
il_6	fixed	NA	pre_postpost	-0.200	0.161	-1.240	11.000	0.241	ns
il_8	fixed	NA	pre_postpost	0.068	0.080	0.842	11.000	0.418	ns
mcp_1	fixed	NA	pre_postpost	0.006	0.047	0.126	11.000	0.902	ns
tn_fa	fixed	NA	pre_postpost	-0.174	0.084	-2.070	11.000	0.063	ns

Significant cytokines

# extract statistically significant cytokines 
cytokines_trtR_log2_psig <- cyto_trtR_log2coefonly %>%
  filter(cyto_trtR_log2coefonly$p.value < 0.05)

print(cytokines_trtR_log2_psig$cytokines)

## [1] "gm_csf"   "il_12p70"

note 11/30/23: Log transforming changes results a bit. IL-5 is no longer significant (pval = 0.053), and now GM-CSF and IL-12p70 are significant. I’ve done nonparametric stats on untransformed data, these 2 cytokines along with IL-5 were significant. The question is, should I be doing log transformation for the mixed linear models? Q-Q plot looks better for IL-5 after log transforming…

Paired t tests

Red

compare_means(log2_cyto_conc_pg_ml ~ pre_post, cytokines_R_log2, method = "t.test", paired = TRUE, group.by = "cytokines")

## # A tibble: 15 × 9
##    cytokines .y.             group1 group2      p p.adj p.format p.signif method
##    <chr>     <chr>           <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr> 
##  1 if_ng     log2_cyto_conc… pre    post   0.147   1    0.147    ns       T-test
##  2 il_1b     log2_cyto_conc… pre    post   0.193   1    0.193    ns       T-test
##  3 il_2      log2_cyto_conc… pre    post   0.350   1    0.350    ns       T-test
##  4 il_6      log2_cyto_conc… pre    post   0.241   1    0.241    ns       T-test
##  5 il_8      log2_cyto_conc… pre    post   0.418   1    0.418    ns       T-test
##  6 il_10     log2_cyto_conc… pre    post   0.847   1    0.847    ns       T-test
##  7 il_12p70  log2_cyto_conc… pre    post   0.0249  0.35 0.025    *        T-test
##  8 mcp_1     log2_cyto_conc… pre    post   0.902   1    0.902    ns       T-test
##  9 tn_fa     log2_cyto_conc… pre    post   0.0628  0.75 0.063    ns       T-test
## 10 il_13     log2_cyto_conc… pre    post   0.540   1    0.540    ns       T-test
## 11 il_5      log2_cyto_conc… pre    post   0.0531  0.69 0.053    ns       T-test
## 12 il_1ra    log2_cyto_conc… pre    post   0.381   1    0.381    ns       T-test
## 13 il_12p40  log2_cyto_conc… pre    post   0.0928  1    0.093    ns       T-test
## 14 gm_csf    log2_cyto_conc… pre    post   0.0124  0.19 0.012    *        T-test
## 15 il_4      log2_cyto_conc… pre    post   0.241   1    0.241    ns       T-test

We get the same result for paired t-test on log transformed data

Immune Cells

Mixed linear modeling

note 11/30/23: if there are any NAs, should I assume below LOD? Or were values not recorded?

# convert immune cell data from wide to long
cells_long <- meta_table %>%
  select(1:24, starts_with("x")) %>%
  pivot_longer(cols = x01_cd45_cd66b_lymph_dc_mono:x41_cd66b_mdsc_grans,
               names_to = "cell_type",
               values_to = "cell_value")

Sequence effect test

# lmer model function
cells_seq_function <- function(cells_long) {
  lmer(cell_value ~ sequence + (1|patient_id), data = cells_long, REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cells_seq_models <- map(split(cells_long,
                              cells_long$cell_type),
                   cells_seq_function)

# create single data frame of results. apply tidy from broom.mixed package to clean up results
cells_seq_results <- map_df(cells_seq_models,
                            tidy,
                            .id = "cell_type")

cells_seq_results_coef <- cells_seq_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Results

kable(cells_seq_results_coef, format = "markdown", digits = 3)

cell_type	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
x01_cd45_cd66b_lymph_dc_mono	fixed	NA	sequenceY_R	-0.024	2.152	-0.011	10.000	0.991	ns
x02_cd45_cd66b_grans	fixed	NA	sequenceY_R	0.445	0.558	0.797	10.000	0.444	ns
x03_cd3_cd45_cd3_t_cells	fixed	NA	sequenceY_R	5.038	7.441	0.677	10.000	0.514	ns
x04_tc_rgd_cd3_ab_t_cells	fixed	NA	sequenceY_R	4.702	7.682	0.612	10.000	0.554	ns
x05_cd4_cd8_cd8_t_cells	fixed	NA	sequenceY_R	3.751	2.477	1.514	10.000	0.161	ns
x06_cd45ro_cd45ra_naive_cd8	fixed	NA	sequenceY_R	1.338	1.832	0.730	10.000	0.482	ns
x07_cd46ro_cd45ra_cm_cd8	fixed	NA	sequenceY_R	0.152	0.210	0.724	10.000	0.485	ns
x08_cd45ro_cd45ra_em_cd8	fixed	NA	sequenceY_R	1.094	1.512	0.723	10.000	0.486	ns
x09_cd45r0_cd45ra_te_cd8	fixed	NA	sequenceY_R	-0.109	0.437	-0.249	10.000	0.809	ns
x10_cd38_hladr_activated_cd8	fixed	NA	sequenceY_R	-0.002	0.026	-0.089	10.000	0.931	ns
x11_cd4_cd8_cd4_t_cells	fixed	NA	sequenceY_R	0.546	7.572	0.072	10.000	0.944	ns
x12_cd45ro_cd45ra_naive_cd4	fixed	NA	sequenceY_R	1.412	4.877	0.290	10.000	0.778	ns
x13_cd45ro_cd45ra_cm_cd4	fixed	NA	sequenceY_R	-0.649	3.423	-0.190	10.000	0.853	ns
x14_cd45ro_cd45ra	fixed	NA	sequenceY_R	0.835	0.737	1.133	10.000	0.284	ns
x15_cd45ro_cd45ra_te_cd4	fixed	NA	sequenceY_R	-0.337	0.565	-0.597	10.001	0.564	ns
x16_cd38_hladr_activated_cd4	fixed	NA	sequenceY_R	0.001	0.087	0.011	10.000	0.991	ns
x17_cd25_cd127_tregs	fixed	NA	sequenceY_R	0.290	0.296	0.978	10.000	0.351	ns
x18_ccr4_cd4_total_ccr4_treg	fixed	NA	sequenceY_R	0.312	0.268	1.163	10.000	0.272	ns
x19_cd45ra_cd45ro_ccr4_treg_naive	fixed	NA	sequenceY_R	0.019	0.011	1.802	10.000	0.102	ns
x20_hladr_total_ccr4_treg_activated	fixed	NA	sequenceY_R	0.008	0.061	0.137	10.000	0.893	ns
x21_cd45ra_cd45ro_ccr4_treg_memory	fixed	NA	sequenceY_R	0.224	0.204	1.097	10.000	0.298	ns
x22_cxcr3_ccr6_th1	fixed	NA	sequenceY_R	0.634	0.708	0.896	10.000	0.391	ns
x23_cxcr3_ccr6_th2	fixed	NA	sequenceY_R	7.379	5.644	1.307	10.000	0.220	ns
x24_cxcr3_ccr6_th17	fixed	NA	sequenceY_R	0.684	2.080	0.329	10.000	0.749	ns
x25_cd19_cd3_b_cells	fixed	NA	sequenceY_R	1.080	2.731	0.395	10.000	0.701	ns
x26_cd27_ig_d_naive_b_cells	fixed	NA	sequenceY_R	1.249	2.409	0.518	10.000	0.615	ns
x27_cd27_ig_d_memory_b_cells	fixed	NA	sequenceY_R	-0.236	0.322	-0.733	10.000	0.480	ns
x28_cd27_ig_d_memory_resting_b_cells	fixed	NA	sequenceY_R	-0.177	0.219	-0.806	10.000	0.439	ns
x30_cd27_cd38_plasmablasts	fixed	NA	sequenceY_R	-0.033	0.021	-1.613	10.000	0.138	ns
x31_cd14_monocytes	fixed	NA	sequenceY_R	-5.481	6.350	-0.863	10.000	0.408	ns
x32_cd16_non_classical_mono	fixed	NA	sequenceY_R	-0.440	0.703	-0.627	10.000	0.545	ns
x33_cd16_classical_mono	fixed	NA	sequenceY_R	-4.477	5.667	-0.790	10.000	0.448	ns
x34_hladr_cd56	fixed	NA	sequenceY_R	-4.413	5.412	-0.815	10.000	0.434	ns
x35_cd16_cd123_cd11c_p_dc	fixed	NA	sequenceY_R	-0.092	0.069	-1.346	10.000	0.208	ns
x36_cd16_cd123_cd11c_m_dc	fixed	NA	sequenceY_R	-3.230	4.565	-0.708	10.000	0.495	ns
x37_cd56_cd161_cd123_nk_cells	fixed	NA	sequenceY_R	-1.323	2.273	-0.582	10.000	0.573	ns
x38_cd16_nk_cells	fixed	NA	sequenceY_R	1.153	2.027	0.569	10.000	0.582	ns
x40_cd14_mdsc_mono	fixed	NA	sequenceY_R	-4.353	3.288	-1.324	10.000	0.215	ns
x41_cd66b_mdsc_grans	fixed	NA	sequenceY_R	-0.035	0.128	-0.274	10.000	0.789	ns

# extract statistically significant cytokines 
cells_seq_psig <- cells_seq_results_coef %>%
  filter(cells_seq_results_coef$p.value < 0.05)

print(cells_seq_psig$cell_type)

## character(0)

*No sequence effects

Intervention comparison

# filter data for post-intervention only
cells_post_long <- cells_long %>%
  filter(intervention != "Baseline") %>%
  filter(pre_post == "post")

str(cells_post_long)

## tibble [936 × 26] (S3: tbl_df/tbl/data.frame)
##  $ patient_id       : Factor w/ 12 levels "6101","6102",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ sex              : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
##  $ age_at_enrollment: num [1:936] 58 58 58 58 58 58 58 58 58 58 ...
##  $ bmi_at_enrollment: num [1:936] 31.1 31.1 31.1 31.1 31.1 ...
##  $ intervention     : Factor w/ 3 levels "Baseline","Red",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ sequence         : Factor w/ 2 levels "R_Y","Y_R": 1 1 1 1 1 1 1 1 1 1 ...
##  $ intervention_week: Factor w/ 5 levels "0","2","6","10",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ pre_post         : Factor w/ 2 levels "pre","post": 2 2 2 2 2 2 2 2 2 2 ...
##  $ period           : Factor w/ 5 levels "B0","B1","B2",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ if_ng            : num [1:936] 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 ...
##  $ il_1b            : num [1:936] 17.9 17.9 17.9 17.9 17.9 ...
##  $ il_2             : num [1:936] 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 ...
##  $ il_6             : num [1:936] 2.79 2.79 2.79 2.79 2.79 2.79 2.79 2.79 2.79 2.79 ...
##  $ il_8             : num [1:936] 3.53 3.53 3.53 3.53 3.53 3.53 3.53 3.53 3.53 3.53 ...
##  $ il_10            : num [1:936] 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 ...
##  $ il_12p70         : num [1:936] 2.79 2.79 2.79 2.79 2.79 2.79 2.79 2.79 2.79 2.79 ...
##  $ mcp_1            : num [1:936] 60.4 60.4 60.4 60.4 60.4 ...
##  $ tn_fa            : num [1:936] 16.8 16.8 16.8 16.8 16.8 ...
##  $ il_13            : num [1:936] 7.47 7.47 7.47 7.47 7.47 7.47 7.47 7.47 7.47 7.47 ...
##  $ il_5             : num [1:936] 3.74 3.74 3.74 3.74 3.74 3.74 3.74 3.74 3.74 3.74 ...
##  $ il_1ra           : num [1:936] 4.29 4.29 4.29 4.29 4.29 4.29 4.29 4.29 4.29 4.29 ...
##  $ il_12p40         : num [1:936] 45.3 45.3 45.3 45.3 45.3 ...
##  $ gm_csf           : num [1:936] 6.46 6.46 6.46 6.46 6.46 6.46 6.46 6.46 6.46 6.46 ...
##  $ il_4             : num [1:936] NA NA NA NA NA NA NA NA NA NA ...
##  $ cell_type        : chr [1:936] "x01_cd45_cd66b_lymph_dc_mono" "x02_cd45_cd66b_grans" "x03_cd3_cd45_cd3_t_cells" "x04_tc_rgd_cd3_ab_t_cells" ...
##  $ cell_value       : num [1:936] 98.61 1.36 36.48 31.92 15.05 ...

# remove x41_cd66b_mdsc_grans because it's causing an error
cells_post_long_edited <- cells_post_long %>%
  filter(cell_type != "x41_cd66b_mdsc_grans")

# treatment effect test lmer model function
cells_trt_function <- function(cells_post_long_edited) {
  lmer(cell_value ~ intervention + (1|patient_id), 
       data = cells_post_long_edited, 
       REML = TRUE)
}

# break data up into subsets based on cell, then apply funtion to each subset
cells_trt_models <- map(split(cells_post_long_edited,
                              cells_post_long_edited$cell_type),
                        cells_trt_function)

# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cells_trt_results <- map_df(cells_trt_models,
                           tidy,
                          .id = "cell_type")

# extract fixed coefficients for sequence only
cells_trt_coefonly <- cells_trt_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

kable(cells_trt_coefonly, format = "markdown", digits = 4)

cell_type	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
x01_cd45_cd66b_lymph_dc_mono	fixed	NA	interventionYellow	0.1437	1.5446	0.0930	11	0.9276	ns
x02_cd45_cd66b_grans	fixed	NA	interventionYellow	-0.3867	0.3053	-1.2669	11	0.2313	ns
x03_cd3_cd45_cd3_t_cells	fixed	NA	interventionYellow	-1.9095	4.0782	-0.4682	11	0.6488	ns
x04_tc_rgd_cd3_ab_t_cells	fixed	NA	interventionYellow	-2.0179	4.0621	-0.4968	11	0.6291	ns
x05_cd4_cd8_cd8_t_cells	fixed	NA	interventionYellow	-0.0913	1.2534	-0.0729	11	0.9432	ns
x06_cd45ro_cd45ra_naive_cd8	fixed	NA	interventionYellow	0.1903	0.6862	0.2773	11	0.7867	ns
x07_cd46ro_cd45ra_cm_cd8	fixed	NA	interventionYellow	-0.0137	0.0717	-0.1914	11	0.8517	ns
x08_cd45ro_cd45ra_em_cd8	fixed	NA	interventionYellow	-0.6942	1.1263	-0.6164	11	0.5502	ns
x09_cd45r0_cd45ra_te_cd8	fixed	NA	interventionYellow	0.3789	0.2313	1.6379	11	0.1297	ns
x10_cd38_hladr_activated_cd8	fixed	NA	interventionYellow	0.0168	0.0119	1.4113	11	0.1858	ns
x11_cd4_cd8_cd4_t_cells	fixed	NA	interventionYellow	-1.6574	2.9437	-0.5630	11	0.5847	ns
x12_cd45ro_cd45ra_naive_cd4	fixed	NA	interventionYellow	-0.6033	1.9806	-0.3046	11	0.7664	ns
x13_cd45ro_cd45ra_cm_cd4	fixed	NA	interventionYellow	-0.1943	0.9081	-0.2140	11	0.8345	ns
x14_cd45ro_cd45ra	fixed	NA	interventionYellow	0.0149	0.2516	0.0594	11	0.9537	ns
x15_cd45ro_cd45ra_te_cd4	fixed	NA	interventionYellow	-1.0865	1.1069	-0.9816	22	0.3370	ns
x16_cd38_hladr_activated_cd4	fixed	NA	interventionYellow	0.0301	0.0376	0.8014	11	0.4398	ns
x17_cd25_cd127_tregs	fixed	NA	interventionYellow	0.1232	0.0949	1.2986	11	0.2207	ns
x18_ccr4_cd4_total_ccr4_treg	fixed	NA	interventionYellow	0.1147	0.0958	1.1970	11	0.2565	ns
x19_cd45ra_cd45ro_ccr4_treg_naive	fixed	NA	interventionYellow	0.0058	0.0069	0.8350	11	0.4215	ns
x20_hladr_total_ccr4_treg_activated	fixed	NA	interventionYellow	0.0527	0.0390	1.3536	11	0.2030	ns
x21_cd45ra_cd45ro_ccr4_treg_memory	fixed	NA	interventionYellow	0.0754	0.0581	1.2975	11	0.2210	ns
x22_cxcr3_ccr6_th1	fixed	NA	interventionYellow	-0.5504	0.6386	-0.8618	11	0.4072	ns
x23_cxcr3_ccr6_th2	fixed	NA	interventionYellow	-0.4578	1.7276	-0.2650	11	0.7959	ns
x24_cxcr3_ccr6_th17	fixed	NA	interventionYellow	1.2526	1.1705	1.0701	11	0.3075	ns
x25_cd19_cd3_b_cells	fixed	NA	interventionYellow	0.1449	0.8671	0.1671	11	0.8703	ns
x26_cd27_ig_d_naive_b_cells	fixed	NA	interventionYellow	0.0560	0.8292	0.0676	11	0.9473	ns
x27_cd27_ig_d_memory_b_cells	fixed	NA	interventionYellow	0.0936	0.2156	0.4343	11	0.6725	ns
x28_cd27_ig_d_memory_resting_b_cells	fixed	NA	interventionYellow	-0.0552	0.0831	-0.6645	11	0.5200	ns
x30_cd27_cd38_plasmablasts	fixed	NA	interventionYellow	0.0185	0.0134	1.3743	11	0.1967	ns
x31_cd14_monocytes	fixed	NA	interventionYellow	0.4091	3.8171	0.1072	11	0.9166	ns
x32_cd16_non_classical_mono	fixed	NA	interventionYellow	-0.3330	0.3102	-1.0735	11	0.3060	ns
x33_cd16_classical_mono	fixed	NA	interventionYellow	0.6321	3.5425	0.1784	11	0.8616	ns
x34_hladr_cd56	fixed	NA	interventionYellow	0.6282	3.2256	0.1947	11	0.8491	ns
x35_cd16_cd123_cd11c_p_dc	fixed	NA	interventionYellow	0.0082	0.0508	0.1606	11	0.8753	ns
x36_cd16_cd123_cd11c_m_dc	fixed	NA	interventionYellow	0.6105	2.9571	0.2064	11	0.8402	ns
x37_cd56_cd161_cd123_nk_cells	fixed	NA	interventionYellow	-4.8544	1.6409	-2.9583	11	0.0130	*
x38_cd16_nk_cells	fixed	NA	interventionYellow	3.3508	0.8126	4.1235	11	0.0017	**
x40_cd14_mdsc_mono	fixed	NA	interventionYellow	-10.9242	3.2788	-3.3317	11	0.0067	**

# extract statistically significant cell type 
cells_trt_psig <- cells_trt_coefonly %>%
  filter(cells_trt_coefonly$p.value < 0.05)


print(cells_trt_psig$cell_type)

## [1] "x37_cd56_cd161_cd123_nk_cells" "x38_cd16_nk_cells"            
## [3] "x40_cd14_mdsc_mono"

*3 significant cell types

Plots

sigcells_trt <- cells_trt_psig$cell_type

lineplots

cells_post_long_edited %>% 
  filter(cell_type %in% sigcells_trt) %>%
  ggplot(aes(x = intervention, y = cell_value, color = patient_id)) +
  geom_line(aes(group = patient_id)) +
  labs(x = "",
       y = "",
       title = "") +
  theme_classic() +
  facet_wrap(~ cell_type, scales = "free_y")

boxplots

cells_post_long_edited %>%
  filter(cell_type %in% sigcells_trt) %>%
  ggplot(aes(x = intervention, y = cell_value, fill = intervention)) +
  geom_boxplot(outlier.shape = NA) +
  scale_fill_manual(values = c("Yellow" = "gold",
                               "Red" = "tomato1")) +
  geom_jitter() +
  labs(x = "Intervention", 
       y = "Cell value",
       title = "Significantly different cell types: post- yellow vs. post- red interventions") +
  facet_wrap(~ cell_type, ncol = 3, nrow = 1, scales = "free_y") +
  theme_bw()

• x38_CD16 NK cells are significantly higher in post- yellow compared to post-red

• x37_cd56_cd161_cd123_nk_cells and x40_cd14_mdsc_mono significantly lower post-yellow intervention compared to post-red

Yellow juice treatment

# subset data into yellow results
cells_Y_long <- cells_long %>%
  filter(intervention == "Yellow")

# remove x41_cd66b_mdsc_grans because it has too many NAs
cells_Y_long_edited <- cells_Y_long %>%
  filter(cell_type != "x41_cd66b_mdsc_grans")

# treatment effect test lmer model function
cells_trtY_function <- function(cells_Y_long_edited) {
  lmer(cell_value ~ pre_post + (1|patient_id), 
       data = cells_Y_long_edited, 
       REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cells_trtY_models <- map(split(cells_Y_long_edited,
                              cells_Y_long_edited$cell_type),
                        cells_trtY_function)

# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cells_trtY_results <- map_df(cells_trtY_models,
                            tidy,
                            .id = "cell_type")

# extract fixed coefficients for sequence only
cells_trtY_coefonly <- cells_trtY_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Results

# extract statistically significant cytokines 
cells_trtY_psig <- cells_trtY_coefonly %>%
  filter(cells_trtY_coefonly$p.value < 0.05)

print(cells_trtY_psig$cell_type)

## [1] "x05_cd4_cd8_cd8_t_cells"     "x08_cd45ro_cd45ra_em_cd8"   
## [3] "x09_cd45r0_cd45ra_te_cd8"    "x15_cd45ro_cd45ra_te_cd4"   
## [5] "x26_cd27_ig_d_naive_b_cells"

Plots

sigcells_Y <- cells_trtY_psig$cell_type

lineplots

cells_Y_long_edited %>% 
  filter(cell_type %in% sigcells_Y) %>%
  ggplot(aes(x = pre_post, y = cell_value, color = patient_id)) +
  geom_line(aes(group = patient_id)) +
  labs(x = "",
       y = "",
       title = "") +
  theme_classic() +
  facet_wrap(~ cell_type, scales = "free_y")

boxplots

cells_Y_long_edited %>% 
  filter(cell_type %in% sigcells_Y) %>%
  ggplot(aes(x = pre_post, y = cell_value, fill = pre_post)) +
  geom_boxplot(outlier.shape = NA) +
  scale_fill_manual(values = c("pre" = "gray",
                               "post" = "gold")) +
  geom_jitter() +
  labs(x = "Pre vs post", 
       y = "Cell value",
       title = "Significantly different cell types: pre vs post yellow") +
  facet_wrap(~ cell_type, ncol = 3, nrow = 2, scales = "free_y") +
  theme_bw()

*5 significant cell types

Red juice treatment

# subset data into red results
cells_R_long <- cells_long %>%
  filter(intervention == "Red")

# remove x41_cd66b_mdsc_grans because it has too many NAs
cells_R_long_edited <- cells_R_long %>%
  filter(cell_type != "x41_cd66b_mdsc_grans")

# treatment effect test lmer model function
cells_trtR_function <- function(cells_R_long_edited) {
  lmer(cell_value ~ pre_post + (1|patient_id), 
       data = cells_R_long_edited, 
       REML = TRUE)
}

# break data up into subsets based on cell, then apply funtion to each subset
cells_trtR_models <- map(split(cells_R_long_edited,
                               cells_R_long_edited$cell_type),
                         cells_trtR_function)

# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cells_trtR_results <- map_df(cells_trtR_models,
                             tidy,
                             .id = "cell_type")

# extract fixed coefficients for sequence only
cells_trtR_coefonly <- cells_trtR_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Results

kable(cells_trtR_coefonly, format = "markdown", digits = 3)

cell_type	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
x01_cd45_cd66b_lymph_dc_mono	fixed	NA	pre_postpost	-0.476	0.978	-0.487	11	0.636	ns
x02_cd45_cd66b_grans	fixed	NA	pre_postpost	0.354	0.228	1.556	11	0.148	ns
x03_cd3_cd45_cd3_t_cells	fixed	NA	pre_postpost	0.294	6.161	0.048	11	0.963	ns
x04_tc_rgd_cd3_ab_t_cells	fixed	NA	pre_postpost	0.086	6.144	0.014	11	0.989	ns
x05_cd4_cd8_cd8_t_cells	fixed	NA	pre_postpost	1.983	1.198	1.655	11	0.126	ns
x06_cd45ro_cd45ra_naive_cd8	fixed	NA	pre_postpost	0.137	0.463	0.296	11	0.773	ns
x07_cd46ro_cd45ra_cm_cd8	fixed	NA	pre_postpost	0.063	0.049	1.296	11	0.221	ns
x08_cd45ro_cd45ra_em_cd8	fixed	NA	pre_postpost	1.645	1.099	1.496	11	0.163	ns
x09_cd45r0_cd45ra_te_cd8	fixed	NA	pre_postpost	0.147	0.191	0.770	11	0.457	ns
x10_cd38_hladr_activated_cd8	fixed	NA	pre_postpost	-0.005	0.014	-0.346	11	0.736	ns
x11_cd4_cd8_cd4_t_cells	fixed	NA	pre_postpost	-2.279	5.045	-0.452	11	0.660	ns
x12_cd45ro_cd45ra_naive_cd4	fixed	NA	pre_postpost	-1.932	2.696	-0.717	11	0.488	ns
x13_cd45ro_cd45ra_cm_cd4	fixed	NA	pre_postpost	-0.585	0.924	-0.633	11	0.540	ns
x14_cd45ro_cd45ra	fixed	NA	pre_postpost	0.374	0.255	1.467	11	0.170	ns
x15_cd45ro_cd45ra_te_cd4	fixed	NA	pre_postpost	1.147	1.105	1.037	22	0.311	ns
x16_cd38_hladr_activated_cd4	fixed	NA	pre_postpost	-0.041	0.031	-1.304	11	0.219	ns
x17_cd25_cd127_tregs	fixed	NA	pre_postpost	-0.063	0.070	-0.903	11	0.386	ns
x18_ccr4_cd4_total_ccr4_treg	fixed	NA	pre_postpost	-0.056	0.058	-0.967	11	0.354	ns
x19_cd45ra_cd45ro_ccr4_treg_naive	fixed	NA	pre_postpost	-0.003	0.004	-0.741	11	0.474	ns
x20_hladr_total_ccr4_treg_activated	fixed	NA	pre_postpost	-0.006	0.025	-0.238	11	0.816	ns
x21_cd45ra_cd45ro_ccr4_treg_memory	fixed	NA	pre_postpost	-0.041	0.047	-0.885	11	0.395	ns
x22_cxcr3_ccr6_th1	fixed	NA	pre_postpost	0.631	0.511	1.234	11	0.243	ns
x23_cxcr3_ccr6_th2	fixed	NA	pre_postpost	-4.914	4.239	-1.159	11	0.271	ns
x24_cxcr3_ccr6_th17	fixed	NA	pre_postpost	-0.076	0.433	-0.175	11	0.864	ns
x25_cd19_cd3_b_cells	fixed	NA	pre_postpost	2.985	0.809	3.690	11	0.004	**
x26_cd27_ig_d_naive_b_cells	fixed	NA	pre_postpost	2.506	0.771	3.250	11	0.008	**
x27_cd27_ig_d_memory_b_cells	fixed	NA	pre_postpost	0.163	0.154	1.055	11	0.314	ns
x28_cd27_ig_d_memory_resting_b_cells	fixed	NA	pre_postpost	0.061	0.028	2.152	11	0.054	ns
x30_cd27_cd38_plasmablasts	fixed	NA	pre_postpost	-0.018	0.018	-1.004	11	0.337	ns
x31_cd14_monocytes	fixed	NA	pre_postpost	-2.571	4.172	-0.616	11	0.550	ns
x32_cd16_non_classical_mono	fixed	NA	pre_postpost	-0.594	0.483	-1.231	11	0.244	ns
x33_cd16_classical_mono	fixed	NA	pre_postpost	-1.970	3.736	-0.527	11	0.608	ns
x34_hladr_cd56	fixed	NA	pre_postpost	-2.146	3.763	-0.570	11	0.580	ns
x35_cd16_cd123_cd11c_p_dc	fixed	NA	pre_postpost	-0.078	0.060	-1.301	11	0.220	ns
x36_cd16_cd123_cd11c_m_dc	fixed	NA	pre_postpost	-1.310	3.138	-0.417	11	0.684	ns
x37_cd56_cd161_cd123_nk_cells	fixed	NA	pre_postpost	-1.292	1.590	-0.812	11	0.434	ns
x38_cd16_nk_cells	fixed	NA	pre_postpost	-0.601	0.521	-1.153	11	0.273	ns
x40_cd14_mdsc_mono	fixed	NA	pre_postpost	-1.314	2.014	-0.652	11	0.528	ns

# extract statistically significant cells 
cells_trtR_psig <- cells_trtR_coefonly %>%
  filter(cells_trtR_coefonly$p.value < 0.05)


print(cells_trtR_psig$cell_type)

## [1] "x25_cd19_cd3_b_cells"        "x26_cd27_ig_d_naive_b_cells"

Plots

sigcells_R <- cells_trtR_psig$cell_type

lineplots

cells_R_long_edited %>% 
  filter(cell_type %in% sigcells_R) %>%
  ggplot(aes(x = pre_post, y = cell_value, color = patient_id)) +
  geom_line(aes(group = patient_id)) +
  labs(x = "",
       y = "",
       title = "") +
  theme_classic() +
  facet_wrap(~ cell_type, scales = "free_y")

boxplots

cells_R_long_edited %>% 
  filter(cell_type %in% sigcells_R) %>%
  ggplot(aes(x = pre_post, y = cell_value, fill = pre_post)) +
  geom_boxplot(outlier.shape = NA) +
  scale_fill_manual(values = c("pre" = "gray",
                               "post" = "tomato")) +
  geom_jitter() +
  labs(x = "Pre vs post", 
       y = "Cell value",
       title = "Significantly different cell types: pre vs post Red") +
  facet_wrap(~ cell_type, ncol = 3, nrow = 2, scales = "free_y") +
  theme_bw()

Naive B-cell pop (#26) significantly increase post-Red intervention, and also increase post-Yellow intervention. This could suggest there is a tomato effect.

Log transformed immune cells

resid normality checks

I need to figure out how to make a function that extracts makes a QQplot for each model. But since IL-5 was significant, let’s look at how the residuals look for this cytokine (for every comparison) before log transforming

#qqplot_fx <- function(cytokine){
  #qqnorm(resid(cyto_seq_models$cytokine))
  #qqline(resid(cyto_seq_models$cytokine))
#}

#qqplot_fx("il_5")

#cyto_seq_results$

seq

qqnorm(resid(cells_seq_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_seq_models$x25_cd19_cd3_b_cells))

Doesn’t look bad

intervention

qqnorm(resid(cells_trt_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_trt_models$x25_cd19_cd3_b_cells))

doesn’t look as bad as sequence

red

qqnorm(resid(cells_trtR_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_trtR_models$x25_cd19_cd3_b_cells))

yellow

qqnorm(resid(cells_trtY_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_trtY_models$x25_cd19_cd3_b_cells))

I will log transform the data and perform mixed linear modeling and paired t-tests for every comparison to see how the data compares.

Wrangle

cells_log2_long <- cells_long %>%
  mutate(log2_cell_value = log2(cell_value)) %>%
  filter(log2_cell_value != -Inf) # remove 0s

Sequence effect test

# lmer model function
cells_log2_seq_function <- function(cells_log2_long) {
  lmer(log2_cell_value ~ sequence + (1|patient_id), data = cells_log2_long, REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cells_log2_seq_models <- map(split(cells_log2_long,
                              cells_log2_long$cell_type),
                   cells_log2_seq_function)

# create single data frame of results. apply tidy from broom.mixed package to clean up results
cells_log2_seq_results <- map_df(cells_log2_seq_models,
                            tidy,
                            .id = "cell_type")

cells_log2_seq_results_coef <- cells_log2_seq_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Results

kable(cells_log2_seq_results_coef, format = "markdown", digits = 3)

cell_type	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
x01_cd45_cd66b_lymph_dc_mono	fixed	NA	sequenceY_R	0.002	0.035	0.064	10.000	0.950	ns
x02_cd45_cd66b_grans	fixed	NA	sequenceY_R	0.084	0.876	0.096	10.000	0.926	ns
x03_cd3_cd45_cd3_t_cells	fixed	NA	sequenceY_R	0.226	0.317	0.712	10.000	0.492	ns
x04_tc_rgd_cd3_ab_t_cells	fixed	NA	sequenceY_R	0.213	0.340	0.625	10.000	0.546	ns
x05_cd4_cd8_cd8_t_cells	fixed	NA	sequenceY_R	0.573	0.367	1.562	10.000	0.149	ns
x06_cd45ro_cd45ra_naive_cd8	fixed	NA	sequenceY_R	0.658	0.858	0.767	10.000	0.461	ns
x07_cd46ro_cd45ra_cm_cd8	fixed	NA	sequenceY_R	0.500	0.516	0.969	10.000	0.355	ns
x08_cd45ro_cd45ra_em_cd8	fixed	NA	sequenceY_R	0.891	0.711	1.254	10.000	0.239	ns
x09_cd45r0_cd45ra_te_cd8	fixed	NA	sequenceY_R	-0.186	0.431	-0.432	10.000	0.675	ns
x10_cd38_hladr_activated_cd8	fixed	NA	sequenceY_R	0.068	0.398	0.172	10.000	0.867	ns
x11_cd4_cd8_cd4_t_cells	fixed	NA	sequenceY_R	0.116	0.507	0.230	10.000	0.823	ns
x12_cd45ro_cd45ra_naive_cd4	fixed	NA	sequenceY_R	0.489	0.725	0.675	10.000	0.515	ns
x13_cd45ro_cd45ra_cm_cd4	fixed	NA	sequenceY_R	0.055	0.515	0.108	10.000	0.916	ns
x14_cd45ro_cd45ra	fixed	NA	sequenceY_R	0.383	0.379	1.013	10.000	0.335	ns
x15_cd45ro_cd45ra_te_cd4	fixed	NA	sequenceY_R	0.964	0.886	1.088	10.000	0.302	ns
x16_cd38_hladr_activated_cd4	fixed	NA	sequenceY_R	0.141	0.408	0.345	10.000	0.737	ns
x17_cd25_cd127_tregs	fixed	NA	sequenceY_R	0.334	0.449	0.745	10.000	0.473	ns
x18_ccr4_cd4_total_ccr4_treg	fixed	NA	sequenceY_R	0.495	0.448	1.105	10.001	0.295	ns
x19_cd45ra_cd45ro_ccr4_treg_naive	fixed	NA	sequenceY_R	0.727	0.688	1.056	10.000	0.316	ns
x20_hladr_total_ccr4_treg_activated	fixed	NA	sequenceY_R	0.198	0.364	0.545	10.000	0.598	ns
x21_cd45ra_cd45ro_ccr4_treg_memory	fixed	NA	sequenceY_R	0.422	0.442	0.956	10.000	0.362	ns
x22_cxcr3_ccr6_th1	fixed	NA	sequenceY_R	1.429	1.201	1.190	10.000	0.262	ns
x23_cxcr3_ccr6_th2	fixed	NA	sequenceY_R	1.660	1.442	1.151	10.000	0.276	ns
x24_cxcr3_ccr6_th17	fixed	NA	sequenceY_R	1.479	1.223	1.209	10.000	0.254	ns
x25_cd19_cd3_b_cells	fixed	NA	sequenceY_R	-0.074	0.456	-0.162	10.000	0.874	ns
x26_cd27_ig_d_naive_b_cells	fixed	NA	sequenceY_R	0.012	0.502	0.023	10.000	0.982	ns
x27_cd27_ig_d_memory_b_cells	fixed	NA	sequenceY_R	-0.461	0.551	-0.837	10.000	0.422	ns
x28_cd27_ig_d_memory_resting_b_cells	fixed	NA	sequenceY_R	-0.400	0.709	-0.564	10.000	0.585	ns
x30_cd27_cd38_plasmablasts	fixed	NA	sequenceY_R	-0.795	0.652	-1.220	9.743	0.251	ns
x31_cd14_monocytes	fixed	NA	sequenceY_R	-0.405	0.448	-0.905	10.000	0.387	ns
x32_cd16_non_classical_mono	fixed	NA	sequenceY_R	-0.531	0.672	-0.790	10.000	0.448	ns
x33_cd16_classical_mono	fixed	NA	sequenceY_R	-0.353	0.447	-0.789	10.000	0.449	ns
x34_hladr_cd56	fixed	NA	sequenceY_R	-0.354	0.442	-0.800	10.000	0.442	ns
x35_cd16_cd123_cd11c_p_dc	fixed	NA	sequenceY_R	-0.757	0.625	-1.211	10.000	0.254	ns
x36_cd16_cd123_cd11c_m_dc	fixed	NA	sequenceY_R	-0.278	0.457	-0.608	10.000	0.557	ns
x37_cd56_cd161_cd123_nk_cells	fixed	NA	sequenceY_R	-0.362	0.644	-0.562	10.000	0.587	ns
x38_cd16_nk_cells	fixed	NA	sequenceY_R	0.213	0.631	0.338	10.000	0.742	ns
x40_cd14_mdsc_mono	fixed	NA	sequenceY_R	-0.711	0.789	-0.901	46.000	0.372	ns
x41_cd66b_mdsc_grans	fixed	NA	sequenceY_R	-0.441	1.071	-0.412	10.000	0.689	ns

# extract statistically significant cytokines 
cells_log2_seq_psig <- cells_log2_seq_results_coef %>%
  filter(cells_log2_seq_results_coef$p.value < 0.05)

print(cells_log2_seq_psig$cell_type)

## character(0)

*No sequence effects

Resid check

# before log transform
qqnorm(resid(cells_seq_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_seq_models$x25_cd19_cd3_b_cells))

# after log transform
qqnorm(resid(cells_log2_seq_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_log2_seq_models$x25_cd19_cd3_b_cells))

Yellow treatment

# subset data into yellow results
cells_log2_Y_long <- cells_log2_long %>%
  filter(intervention == "Yellow")

# lmer model function
cells_log2_Y_function <- function(cells_log2_Y_long) {
  lmer(log2_cell_value ~ sequence + (1|patient_id), data = cells_log2_Y_long, REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cells_log2_Y_models <- map(split(cells_log2_Y_long,
                              cells_log2_Y_long$cell_type),
                   cells_log2_Y_function)

# create single data frame of results. apply tidy from broom.mixed package to clean up results
cells_log2_Y_results <- map_df(cells_log2_Y_models,
                            tidy,
                            .id = "cell_type")

cells_log2_Y_results_coef <- cells_log2_Y_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Results

kable(cells_log2_Y_results_coef, format = "markdown", digits = 3)

cell_type	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
x01_cd45_cd66b_lymph_dc_mono	fixed	NA	sequenceY_R	0.053	0.049	1.086	10	0.303	ns
x02_cd45_cd66b_grans	fixed	NA	sequenceY_R	-0.359	0.994	-0.362	10	0.725	ns
x03_cd3_cd45_cd3_t_cells	fixed	NA	sequenceY_R	-0.072	0.320	-0.226	10	0.826	ns
x04_tc_rgd_cd3_ab_t_cells	fixed	NA	sequenceY_R	-0.109	0.349	-0.311	10	0.762	ns
x05_cd4_cd8_cd8_t_cells	fixed	NA	sequenceY_R	0.488	0.392	1.245	10	0.242	ns
x06_cd45ro_cd45ra_naive_cd8	fixed	NA	sequenceY_R	0.001	0.946	0.001	10	1.000	ns
x07_cd46ro_cd45ra_cm_cd8	fixed	NA	sequenceY_R	0.160	0.533	0.299	10	0.771	ns
x08_cd45ro_cd45ra_em_cd8	fixed	NA	sequenceY_R	1.258	0.713	1.764	10	0.108	ns
x09_cd45r0_cd45ra_te_cd8	fixed	NA	sequenceY_R	-0.085	0.563	-0.151	10	0.883	ns
x10_cd38_hladr_activated_cd8	fixed	NA	sequenceY_R	0.366	0.510	0.717	10	0.490	ns
x11_cd4_cd8_cd4_t_cells	fixed	NA	sequenceY_R	-0.325	0.518	-0.627	10	0.544	ns
x12_cd45ro_cd45ra_naive_cd4	fixed	NA	sequenceY_R	-0.163	0.776	-0.210	10	0.838	ns
x13_cd45ro_cd45ra_cm_cd4	fixed	NA	sequenceY_R	-0.347	0.528	-0.658	10	0.525	ns
x14_cd45ro_cd45ra	fixed	NA	sequenceY_R	0.469	0.424	1.107	10	0.294	ns
x15_cd45ro_cd45ra_te_cd4	fixed	NA	sequenceY_R	1.360	0.860	1.582	10	0.145	ns
x16_cd38_hladr_activated_cd4	fixed	NA	sequenceY_R	0.204	0.515	0.396	10	0.700	ns
x17_cd25_cd127_tregs	fixed	NA	sequenceY_R	-0.016	0.504	-0.033	10	0.975	ns
x18_ccr4_cd4_total_ccr4_treg	fixed	NA	sequenceY_R	0.248	0.519	0.478	10	0.643	ns
x19_cd45ra_cd45ro_ccr4_treg_naive	fixed	NA	sequenceY_R	0.540	0.754	0.716	10	0.490	ns
x20_hladr_total_ccr4_treg_activated	fixed	NA	sequenceY_R	0.032	0.450	0.070	10	0.945	ns
x21_cd45ra_cd45ro_ccr4_treg_memory	fixed	NA	sequenceY_R	0.140	0.497	0.281	10	0.784	ns
x22_cxcr3_ccr6_th1	fixed	NA	sequenceY_R	2.235	1.263	1.769	10	0.107	ns
x23_cxcr3_ccr6_th2	fixed	NA	sequenceY_R	1.188	1.581	0.752	10	0.470	ns
x24_cxcr3_ccr6_th17	fixed	NA	sequenceY_R	0.726	1.201	0.604	10	0.559	ns
x25_cd19_cd3_b_cells	fixed	NA	sequenceY_R	-0.028	0.425	-0.066	10	0.949	ns
x26_cd27_ig_d_naive_b_cells	fixed	NA	sequenceY_R	0.000	0.471	0.000	10	1.000	ns
x27_cd27_ig_d_memory_b_cells	fixed	NA	sequenceY_R	-0.228	0.570	-0.400	10	0.698	ns
x28_cd27_ig_d_memory_resting_b_cells	fixed	NA	sequenceY_R	-0.040	0.706	-0.057	10	0.956	ns
x30_cd27_cd38_plasmablasts	fixed	NA	sequenceY_R	-0.650	0.737	-0.881	10	0.399	ns
x31_cd14_monocytes	fixed	NA	sequenceY_R	0.164	0.489	0.334	10	0.745	ns
x32_cd16_non_classical_mono	fixed	NA	sequenceY_R	0.448	0.711	0.630	10	0.543	ns
x33_cd16_classical_mono	fixed	NA	sequenceY_R	0.217	0.498	0.436	10	0.672	ns
x34_hladr_cd56	fixed	NA	sequenceY_R	0.144	0.482	0.298	10	0.772	ns
x35_cd16_cd123_cd11c_p_dc	fixed	NA	sequenceY_R	-0.232	0.599	-0.387	10	0.707	ns
x36_cd16_cd123_cd11c_m_dc	fixed	NA	sequenceY_R	0.258	0.507	0.509	10	0.622	ns
x37_cd56_cd161_cd123_nk_cells	fixed	NA	sequenceY_R	-0.381	0.724	-0.526	10	0.610	ns
x38_cd16_nk_cells	fixed	NA	sequenceY_R	0.454	0.679	0.670	10	0.518	ns
x40_cd14_mdsc_mono	fixed	NA	sequenceY_R	-0.114	0.923	-0.123	10	0.904	ns

# extract statistically significant cytokines 
cells_log2_Y_psig <- cells_log2_Y_results_coef %>%
  filter(cells_log2_Y_results_coef$p.value < 0.05)

print(cells_log2_Y_psig$cell_type)

## character(0)

Now there arent any significant cells for yellow after log transforming

Resid check

# before log transform
qqnorm(resid(cells_trtY_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_trtY_models$x25_cd19_cd3_b_cells))

# after log transform
qqnorm(resid(cells_log2_Y_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_log2_Y_models$x25_cd19_cd3_b_cells))

Red treatment

# subset data into yellow results
cells_log2_R_long <- cells_log2_long %>%
  filter(intervention == "Red")

# lmer model function
cells_log2_R_function <- function(cells_log2_R_long) {
  lmer(log2_cell_value ~ sequence + (1|patient_id), data = cells_log2_R_long, REML = TRUE)
}

# break data up into subsets based on cytokine, then apply funtion to each subset
cells_log2_R_models <- map(split(cells_log2_R_long,
                              cells_log2_R_long$cell_type),
                   cells_log2_R_function)

# create single data frame of results. apply tidy from broom.mixed package to clean up results
cells_log2_R_results <- map_df(cells_log2_R_models,
                            tidy,
                            .id = "cell_type")

cells_log2_R_results_coef <- cells_log2_R_results %>%
  filter(effect == "fixed") %>%
  filter(term != "(Intercept)") %>%
  add_significance(p.col = "p.value",
                   output.col = "p.sig",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns")
                   )

Results

kable(cells_log2_R_results_coef, format = "markdown", digits = 3)

cell_type	effect	group	term	estimate	std.error	statistic	df	p.value	p.sig
x01_cd45_cd66b_lymph_dc_mono	fixed	NA	sequenceY_R	-0.049	0.032	-1.539	10.00	0.155	ns
x02_cd45_cd66b_grans	fixed	NA	sequenceY_R	0.527	0.808	0.652	10.00	0.529	ns
x03_cd3_cd45_cd3_t_cells	fixed	NA	sequenceY_R	0.524	0.345	1.519	10.00	0.160	ns
x04_tc_rgd_cd3_ab_t_cells	fixed	NA	sequenceY_R	0.534	0.365	1.463	10.00	0.174	ns
x05_cd4_cd8_cd8_t_cells	fixed	NA	sequenceY_R	0.658	0.367	1.792	10.00	0.103	ns
x06_cd45ro_cd45ra_naive_cd8	fixed	NA	sequenceY_R	1.315	0.788	1.670	10.00	0.126	ns
x07_cd46ro_cd45ra_cm_cd8	fixed	NA	sequenceY_R	0.841	0.543	1.550	10.00	0.152	ns
x08_cd45ro_cd45ra_em_cd8	fixed	NA	sequenceY_R	0.525	0.743	0.707	10.00	0.496	ns
x09_cd45r0_cd45ra_te_cd8	fixed	NA	sequenceY_R	-0.287	0.318	-0.903	10.00	0.388	ns
x10_cd38_hladr_activated_cd8	fixed	NA	sequenceY_R	-0.229	0.346	-0.662	10.00	0.523	ns
x11_cd4_cd8_cd4_t_cells	fixed	NA	sequenceY_R	0.558	0.527	1.058	10.00	0.315	ns
x12_cd45ro_cd45ra_naive_cd4	fixed	NA	sequenceY_R	1.141	0.723	1.578	10.00	0.146	ns
x13_cd45ro_cd45ra_cm_cd4	fixed	NA	sequenceY_R	0.458	0.533	0.860	10.00	0.410	ns
x14_cd45ro_cd45ra	fixed	NA	sequenceY_R	0.298	0.369	0.806	10.00	0.439	ns
x15_cd45ro_cd45ra_te_cd4	fixed	NA	sequenceY_R	0.568	1.041	0.546	10.00	0.597	ns
x16_cd38_hladr_activated_cd4	fixed	NA	sequenceY_R	0.078	0.376	0.207	10.00	0.840	ns
x17_cd25_cd127_tregs	fixed	NA	sequenceY_R	0.685	0.418	1.638	10.00	0.132	ns
x18_ccr4_cd4_total_ccr4_treg	fixed	NA	sequenceY_R	0.743	0.410	1.813	10.00	0.100	ns
x19_cd45ra_cd45ro_ccr4_treg_naive	fixed	NA	sequenceY_R	0.914	0.860	1.062	10.00	0.313	ns
x20_hladr_total_ccr4_treg_activated	fixed	NA	sequenceY_R	0.365	0.324	1.126	10.00	0.286	ns
x21_cd45ra_cd45ro_ccr4_treg_memory	fixed	NA	sequenceY_R	0.705	0.409	1.724	10.00	0.115	ns
x22_cxcr3_ccr6_th1	fixed	NA	sequenceY_R	0.623	1.191	0.523	10.00	0.613	ns
x23_cxcr3_ccr6_th2	fixed	NA	sequenceY_R	2.133	1.316	1.620	10.00	0.136	ns
x24_cxcr3_ccr6_th17	fixed	NA	sequenceY_R	2.233	1.285	1.738	10.00	0.113	ns
x25_cd19_cd3_b_cells	fixed	NA	sequenceY_R	-0.120	0.506	-0.237	10.00	0.818	ns
x26_cd27_ig_d_naive_b_cells	fixed	NA	sequenceY_R	0.023	0.565	0.041	10.00	0.968	ns
x27_cd27_ig_d_memory_b_cells	fixed	NA	sequenceY_R	-0.695	0.554	-1.253	10.00	0.239	ns
x28_cd27_ig_d_memory_resting_b_cells	fixed	NA	sequenceY_R	-0.759	0.741	-1.025	10.00	0.330	ns
x30_cd27_cd38_plasmablasts	fixed	NA	sequenceY_R	-0.797	0.531	-1.503	8.25	0.170	ns
x31_cd14_monocytes	fixed	NA	sequenceY_R	-0.974	0.465	-2.096	10.00	0.063	ns
x32_cd16_non_classical_mono	fixed	NA	sequenceY_R	-1.510	0.733	-2.061	10.00	0.066	ns
x33_cd16_classical_mono	fixed	NA	sequenceY_R	-0.922	0.458	-2.013	10.00	0.072	ns
x34_hladr_cd56	fixed	NA	sequenceY_R	-0.851	0.459	-1.854	10.00	0.093	ns
x35_cd16_cd123_cd11c_p_dc	fixed	NA	sequenceY_R	-1.282	0.684	-1.874	10.00	0.090	ns
x36_cd16_cd123_cd11c_m_dc	fixed	NA	sequenceY_R	-0.813	0.463	-1.757	10.00	0.110	ns
x37_cd56_cd161_cd123_nk_cells	fixed	NA	sequenceY_R	-0.343	0.671	-0.512	10.00	0.620	ns
x38_cd16_nk_cells	fixed	NA	sequenceY_R	-0.028	0.629	-0.044	10.00	0.965	ns
x40_cd14_mdsc_mono	fixed	NA	sequenceY_R	-1.309	0.717	-1.825	10.00	0.098	ns
x41_cd66b_mdsc_grans	fixed	NA	sequenceY_R	-0.441	1.071	-0.412	10.00	0.689	ns

# extract statistically significant cytokines 
cells_log2_R_psig <- cells_log2_R_results_coef %>%
  filter(cells_log2_R_results_coef$p.value < 0.05)

print(cells_log2_R_psig$cell_type)

## character(0)

Now there arent any significant cells for red after log transforming

Resid check

# before log transform
qqnorm(resid(cells_trtR_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_trtR_models$x25_cd19_cd3_b_cells))

# after log transform
qqnorm(resid(cells_log2_R_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_log2_R_models$x25_cd19_cd3_b_cells))

Looks even worse after log transform